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beneath  an  embankment 
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r-         -    dimerisiori  used  to  calculate  stress  charige 

beneath  an  embankment 
r^        -    dimension  used  to  calculate  stress  change 

beneath  an  embankment 

R         -    radius 

R  -    strength 

R  .        -    minimum  value  of  strength 
mm  ^ 

S  -    load 

S         -    settlement 

S  -    forward  speed 

-  the  actual  consolidation  settlement  of  a 
compressible  layer 

-  calculated  value  of  consolidation  settlement 
of  a  compressible  layer  assuming  excess 
pore  pressures  in  the  field  equal  to 

the  pressure  developed  in  laboratory 

S  -    maximum  value  of  load 

max 

t  -    time 

T  -    chord  length  of  segments  circumscribing  a 

circle 
T  -    time  to  travel  the  distance  L 

u  -    excess  pore  pressure 

u.   .       -    excess  pore  pressure  at  the 

.th       ...     ..    .th 

1    X  position,  the  J    z  position 

and  the  k    time  step 

UU         -    unconsol idated-undrained  triaxial  test 

U%         -    percent  consolidation 
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the  percent  volume  change  due  to  saturation 

at  the  depth,  z 

variable  used  to  compute  the  friction-circle 

factor  of  safety 

parameter  used  to  defirie  the  beta  distribut  ion 

mold  volume 

variance  of  x 

frequency  of  vibration  of  a  vibratory  compactor 

(ja t er  c  on t  en t 

compactor  weight 

the  coordinate  value  on  the  x  axis 

parameter  u£ed  to  define  the  beta  distrbution 

>'.    coordinate  of  the  intersection  of  the  first 

chord  and  its  perpendicular  bisector 

X  coordinate  of  the  intersection  of  the  second 

chord  and  its  perpendicular  bisector 

X  coordinate  of  the  center  of  the  circle 

X  coordinate  of  the  first  point  on  the  first 

chordofacircle 

X    coordinate    of    the    second   point    ori   the    first 

chord  of  a  circle 

X  coordinate  of  the  second  point  on  the  second 

chord  of  a  circle 

the  coordinate  value  on  the  y  axis 

vertical  distance  to  bottom  of  a  slice 

from  the  moment  center 
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y         -    y  coordiT^«te  of  th«  intersection  of  the  first 

chord  and  its  perpendicular  bisector 

u  _       -    y  coordinate  of  the  intersection  of  the  second 
^c2  ' 

chord  and  its  perpendicular  bisector 
y         -    y  coordinate  of  the  center  of  the  circle 
y         -    y  coordinte  of  the  first  point  on  the  first 

chord  of  a  circle 
y_         ~    y  coordinate  of  the  second  point  on  the  first 

chord  of  a  circle 
y_        -    y  coordinate  of  the  second  point  on  the  second 

chord  of  a  circle 
Y         -   arbitrary  function 
Z         -    net  uork  per  vibratory  cycle  of  a  vibratory 

compactor 
a         -    angle  the  bottom  of  a  slice  makes  with 

the  horizontal 
oc         -    parameter  used  to  define  the  beta  distribution 
a  -    coefficient  used  in  equation  4.19 

a  -    coefficient  used  in  equation  4.19 

B  -        sideslope  angle 

B  -        value  of  sideslope  at  which  a  slope 

reaches  limit  equilibrium 

^C        -    cohesion  force  acting  on  the  bottom  of  a  slice 
a 

Ae        -    the  difference  between  the  initial 
o 

void  ratio  in  the  soil  sample 
and  another  soil  stratum 

AN         -    normal  force  acting  on  the  bottom  of  a  slice 
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&Q        -    surcharge  load  on  a  slice 

AS         -    resisting  force  acting  on  the 

bottom  of  a  slice 
6t        -   time  increment 
6M  ~         slice  weight 

AU         -    uiater  force  acting  beneath 

a  slice 
AU-       -    water  force  acting  on  top  of 

a  slice 
Ax         -    increment  in  the  x  direction 
Az        -    increment  in  the  z  direction 
Az.        -    vertical  grid  spacing  above  a 

layer  interface 
Az-       -    vertical  grid  spacing  beneath  a 

layer  interface 
Aa         -    stress  change 
Aa        -    change  in  vertical  stress 
Act        -    change  in  normal  stress  acting 

in  the  x  direction 
Aa        -    change  in  normal  stress  acting 

in  the  z  direction 
Act-        -   change  in  major  principal  stress 
ACTg       -    change  in  minor  principal  stress 
At        -    change  in  shear  stress  acting  in  the 

z  direction  on  the  plane  perpendicular 

to  the  X  axis 
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&6  -  deflection  angle  of  segments  circtmscribing 
a  circle 

V  -  unit  weight 

V  .  -  drg  unit  mei^t 

V  .  -  maximum  dry  unit  weight 

V_  -  moist  unit  qieight 
m 

V  -  unit  weight  of  water 
r  -  gamma  function 

X  -  nond i mens i ona 1 i zed  slope  stability  parameter 

\k  -  consolidation  settlement  correction  factor 

\i.  -  coefficient  of  variation  of  the 
cohesion  intercept 

V  -  Poisson's  ratio 
•♦  -  friction  angle 

^  -  maximum  available  friction  angle 

♦  -  friction  angle  required  for  equilibrium 

a^^4.  -  octahedral  normal  stress 

OCX 

c''  -  effective  preconsol  idat  ion  pr^sstjrc 

^samole  ~  **^®  vertical  pressure  at  the  depth 

of  the  soil  sample 

cf'^  -  effective  vertical  pressure 

o'"  -  the  effective  overburden  pres5»ure 

in  a  soil  stratum 

o  -  yield  strength 

*r.  -  major  principal  stress 

^2  ~  intermediate  principal  stress 


a  -        minor  principal  stress 
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E  -    summation 

T   ^       -    octahedral  shear  stress 
oct 


HIGHLIGHT  SUMMARY 

This  study  illustrates  how  the  results  of  the  compacted 
clay  investigation  can  be  best  used  in  the  design  and  sta- 
bility analysis  of  compacted  clay  embankments.   It  also  com- 
pletes the  analysis  package  by  supplying  computer  programs 
for  the  calculation  of  embankment  settlement. 

The  alternatives  of  specifying  compaction  procedures  or 
compaction  results  are    compared  and  an  hybrid  approach  of 
specifying  compaction  that  integrates  the  advantages  of 
these  two  approaches  is  introduced.    In  the  hybrid  approach 
of  compaction  specification/  the  compactive  effort  is  speci- 
fied so  that  the  corresponding  optimum  moisture  content  is 
equal  to  the  expected  compaction  water  content. 

Embankment  side  slope  designs  are  illustrated  for  short 
and  long  term  conditions  using  laboratory  compacted  shear 
strength  data.   In  these  examples  the  embankment  material  is 
assumed  to  be  compacted  St.  Croix  clay  and  the  strength 
parameters  for  short  and  long  term  conditions  are  obtained 
from  the  reports  by  Weitzel  and  Lovell  (1979)  and  Johnson 
and  Lovell  {1979)i  respectively.   Several  improvements  made 
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to  the  STABL  program  during  the  course  of  this  study  are 
presented : 

-  The  Simplified  Bishop  factor  of  safety  option  was 
receded  to  correct  various  difficulties  discovered  by 
STABL  users. 

-  Recommendations  are   made  to  avoid  common  errors  in  the 
use  of  STABL  (surfaces  u/ith  unacceptable  shapes,  direc- 
tion limits  on  benched  slopes,  direction  of  surface 
generation,  etc.  ). 

-  A  methodology  was  developed  to  adjust  the  Simplified 
Janbu  factor  of  safety  (used  in  several  STABL  options) 
to  more  familiar  definitions  of  the  factor  of  safety. 

Geometric  and  probabilistic  interpretation  of  the  factor  of 
safety  are  introduced  to  demonstrate  their  usefulness  as 
alternatives  and/or  supplements  to  the  conventional  strength 
factor  of  safety.   The  probabilistic  approach  takes  into 
account  the  variability  of  material  parameters  and  provide 
the  reliability  of  the  slope  corresponding  to  the  computed 
factor  of  safety. 

Computer  programs  are    developed  to  compute  the  magni- 
tude of  settlement  that  occur  within  the  embankment  itself 
as  well  as  of  the  consolidation  settlement  of  fine  grained 
soil  layers  belou;  the  embankment.   Programs  are   also  pro- 
vided to  estimate  the  time-rate  of  consolidation  settlement 
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of  the  compressible  soil  layers  beneath  an  embankment. 
These  programs;  in  conjunction  with  STABL,  form  an  analysis 
package  for  the  design  of  embankments.   User's  manuals  and 
listings  are    given  for  all  these  computer  programs. 


I  -  INTRODUCTION 

This  report  is  a  part  of  Purdue's  embankment  perfor- 
mance project  which  was  initiated  to  improve  the  ability  of 
highujay  engineers  to  design  highway  embankments. 

This  project  is  composed  of  two  separate  objectives. 
One  objective  is  the  investigation  of  the  strength  and 
compressibility  parameters  of  compacted  St.  Croix  clay;  a 
highly  plastic  residual  clay  of  Southern  Indiana.  The  second 
objective  of  this  project  entails  the  development  of  an 
analysis  package  to  predict  the  settlement  and  stability 
performance  of  embankments. 

Compacted  Clau  Investigation 

The  study  of  the  parameters  defining  the  behavior  of 
compacted  St.  Croix  clay  was  performed  in  two  phases.   The 
first  phase  consisted  of  tests  performed  on  clay  samples 
that  were  prepared  to  simulate  standard  laboratory  compac- 
tion specifications.   This  testing  program  included  the  fol- 
lowing separate  investigations: 

1)  Compressibility  and  prestress  behavior  of  St. 
Croix  clay  (DiBernardo,  1979). 


2)  Unconsol idated-undrained  <UU)  shear  strength  behav- 
ior of  St  Croix  clay  (Weitzel,  1979). 

3)  Consol idated-undrained  (CU)  shear  strength  of  St. 
Croix  clay  (Johnson*  1979). 

Parameters  defining  the  behavior  of  field  compacted 
clay  are  different  than  those  obtained  with  laboratory  test- 
ing because  of  the  difference  in  compactive  effort  between 
field  compactors  and  laboratory  compaction  tests.   There- 
fore; in  the  second  phase*  the  tests  of  the  first  phase  were 
repeated  on  samples  compacted  in  the  field  with  two  dif- 
ferent rollers  (Lin.  1981  and  Liang,  19S2>.  The  swell  pres- 
sure of  both  laboratory  and  field  compacted  clay  was  studied 
by  Terdich  (1981).   These  studies  provide  a  unique  look  at 
the  correlation  between  the  behavior  of  field  and  laboratory 
c  ompac  ted  clays. 

Analusis  Package 

Purdue  University  has  a  long  standing  interest  in  the 
development  of  user-oriented  slope  stability  computer  pro- 
grams   One  of  the  first  developments  was  the  SLOPE  program 
package  (Carter.  1971)  consisting  of  four  separate  programs. 
The  subsequent  development  was  the  STABL  program  (Siegel. 
1975).   This  program  can  evaluate  the  factor  of  safety  of 
slopes  of  almost  any  description  and  shape.   Boundary  sur- 
charge loads  and  pseudo-static  earthquake  forces  may  also  be 


included.   The  most  significant  feature  of  STAPL  is  its 

ability  to  automatically  create  randomly  generated  surfaces 

to  help  the  user  search  for  the  minimum  value  of  the  factor 
of  safety. 

STABL  was  further  developed  by  Doutrup  (1977).   Her 
improvements  included: 

1)  The  addition  of  specialized  search  routines  for 
simulation  of  block  shaped  failure  surfaces  with 
active  and  passive  wedges. 

2)  The  inclusion  of  the  Simplified  Bishop  factor  of 
safety. 

3)  Changes  in  the  input  of  pore  water  pressure 
that  allow  the  simulation  of  artesian  pressure. 

S1ABL  is  used  on  a  regular  basis  by  the  Indiana  Department 
of  Highways  for  routine  evaluation  of  slope  stability. 

The  most  recent  of  Purdue  contributions  to  the  field  of 
slope  stability  is  the  development  of  a  factor  of  safety 
that  takes  into  account  the  three-dimensional  nature  of 
limit  equilibrium  surfaces  (Cheni  1981).   Currentlyj  this 
method  is  programmed  only  for  simple  slope  shapes  and 
failure  surfaces. 


Report  Organization 

The  purpose  of  this  report  is  twofold    Firsti  it  com- 
pletes the  analysis  package  by  supplying  computer  programs 
for  the  calculation  of  embankment  settlement.   Secondi  it 
illustrates  hou/  the  results  of  the  compacted  clay  investiga- 
tion may  be  best  used  in  the  design  of  compacted  clay 
embankments. 

Chapter  II  provides  an  overview  of  compaction  specifi- 
cation.  The  alternatives  of  specifying  compaction  pro- 
cedures or  compaction  results  are    compared.   A  hybrid 
approach  of  specifying  compaction  that  integrates  the  advan- 
tages of  these  two  approaches  is  introduced. 

Chapter  III  covers  a  variety  of  topics  pertinent  to  the 
subject  of  slope  stability  including: 

1)  Use  of  the  STABL  program  for  slope  stability  analy- 
sis.  The  discussion  includes  recent  corrections 
made  to  the  Simplified  Bishop  option. 

2)  Strength  parameters  of  compacted  clays  to  be  used 
when  assessing  the  factor  of  safety  and  stability  of 
compacted  clay  embankments 

3)  Discussion  of  geometric  and  probabilistic 
interpretations  of  the  factor  of  safety. 


Chapter  IV  deals  uith  the  settlements  caused  by  the 
construction  of  an  embankment.   Usei — oriented  computer  pro- 
grams to  facilitate  the  computation  of  the  magnitude  and 
time-rate  of  consolidation  settlement  are    included. 


Chapter  V  presents  conclusions  of  the  uork  that  was 
done  and  suggestions  for  further  research. 


II  -  COMPACTION  SPECIFICATION 

Compaction  is  the  densi f i cat i on  of  soil  by  the  applica- 
tion of  mechanical  energy.   Densi f icat i on  is  achieved  by 
reduction  of  the  size  and  number  of  air  voids  in  the  soil. 
As  the  volume  of  the  voids  is  reduced,  the  shear  strength  of 
the  soil  is  increased  and  its  permeability  is  decreased. 
The  increase  in  shear  strength  and  the  reduction  in  permea- 
bility are    two  factors  that  make  compaction  a  good  technique 
for  constructing  highway  embankments  and  dams. 

Specifying  an  adequate  level  of  coiTipaction  and  range  of 
water  content  indirectly  assures  that  the  soil  will  have  a 
relatively  high  shear  strength  and  a  low  compressibility. 
This  is  fortunate  because  these  soil  properties  are    diffi- 
cult and  time  consuming  to  measure  on  a  routine  basis. 

It  is  important  to  be  able  to  quantify  the  level  of 
compaction  to  determine  if  the  compaction  is  adequate.    In 
his  pioneering  work.  Proctor  showed  that  the  level  of  com- 
paction depends  on  the  compactive  effort/  E,  i.  e.  ,     the 
amount  of  mechanical  energy  imparted  into  a  unit  volume  of 


soil  (Proctor.  1933).   For  a  given  compactive  effort,  the 
dry  density  ,  V  ,  that  can  be  achieved  varies  with  the  uiater 
content  of  the  soil  (Figure  2.1).   There  is  a  value  of  the 
biater  content  called  the  optimum  moisture  content  .  OMC.  at 
which  the  dry  density  has  a  maximum  value.   The  soil  shear 
strength  is  maximized  and  the  soil  permeability  in  service 
is  minimized  at  or  near  this  water  content.  The  OMC  is  a 
function  of  the  compactive  effort  that  varies  along  a  "line 
of  optimums"  (Figure  2.2).  In  general,  the  dry  density  will 
increase  and  the  OMC  will  decrease  as  the  compactive  effort 
increases. 

The  concept  of  compactive  effort  provides  the  basis  for 
evaluating  the  level  of  compaction.  At  a  specified  value  of 
compactive  effort,  there  is  a  maximum  achievable  dry  density 

,  v.     ,  corresponding  to  the  OMC  which  can  be  achieved. 

d  max         ^  ^ 

At  water  contents  other  than  the  OMC,  V  .  that  is  achieved 

d 

will  be  less  than  V_.        Therefore,  it  is  convenient  to 
d  max 

define  the  "percent  compaction"  as  the  ratio  of  V  .  to 
d  max 

Specification  of  Procedure 

The  basic  philosophy  of  specifying  the  compactive  pro- 
cedure is  that  if  a  certain  procedure  is  followed,  the  com- 
paction is  assumed  to  be  satisfactory.  Typically.  this  is 
achieved  by  construction  of  a  test  pad.  Each  time  the  com- 
pactor passes  over  a  soil  lift  in  the  test  pad.  the  soil 
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FIGURE     2. 1      MOISTURE-DENSITY  RELATIONSHIP 
OF  A  COMPACTED  CLAY 
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FIGURE  2.2 


VARIATION  OF  THE  MOISTURE  -  DENSITY  RELATIONSHIP  OF 
A  COMPACTED  CLAY  VS.  COMPACTIVE  EFFORT 


10 

density  is  recorded  using  either  a  sand  cone  test  or  a 
nuclear  density  device.   The  measured  density  is  then  plot- 
ted versus  the  number  of  passes  of  the  conpactor.   The 
result  is  a    "density  grouith  curve"  (Figure  2.3).   The  den- 
sity growth  curve  indicates  the  effectiveness  of  each  pass 
of  the  compactor.   The  increase  in  density  diminishes  with 
each  pass  of  the  compactor  until  a  pass  generates  a  negligi- 
ble density  increase.   Beyond  this  point/  higher  densities 
can  only  be  obtained  by  using  machines  that  impart  more  com- 
pactive  effort.   Therefore/  although  the  number  of  passes 
should  be  specified  to  provide  a  large  fraction  (perhaps 
niorc*  thai)  V3'/.  >  of  the  achievable  density/  the  specified 
number  of  passes  should  be  less  than  the  value  at  which  each 
subsequent  pass  creates  only  a  negligible  increase  in  den- 
sity. 

It  is  recommended  that  the  density  growth  curve  be 
developed  for  various  lift  thicknesses/  various  rates  of 
advance  of  the  compactor/  and  different  compactors.   This 
makes  it  possible  to  determine  the  best  compactor  and  the 
optimum  mode  of  operation  for  a  given  job.  The  performance 
study  done  at  Purdue  (Terdich/  1961)  is  a  good  example  of 
the  use  of  the  density  growth  curve.   This  study  showed  that 
a  Caterpillar  Model  825  tamping  compactor  is  significantly 
more  effective  than  a  RayGo  Rascal  Model  420C  Vibratory  com- 
pactor to  compact  St.  Croix  clay.   Advantages  of  specifying 
the  procedure  include: 
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1)  Only  limited  testing  is  required  to  provide  quality 
control  of  the  compaction. 

2)  The  contractor  has  the  assurance  that  if  he  operates 
in  accordance  with  the  specified  procedure*  he 

will  obtain  the  necessary  compaction.  This  helps 
reduce  the  adversarial  aspect  of  the  relationship 
between  the  contractor  and  the  engineer. 


Specification  of  Results 

On  jobs  where  the  compacted  soil  is  expected  to  be 
quite  variable,  specifying  the  procedure  may  not  be  practi- 
cal because  the  level  of  compaction  will  change  in  an  unk- 
nown fashion  from  soil  to  soili  even  if  identical  procedures 
are  used.   When  this  occurs*  the  results  of  the  compactive 
bperation  must  be  specified. 

Specification  of  results  entails  two  stages.   The  first 
stage  involves  the  development  of  a  moisture-density  rela- 
tionship from  compaction  tests  run  on  soil  samples  compacted 
in  the  laboratory.   The  laboratory  compaction  should  produce 
similar  moisture-density  results  as  those  produced  by  the 
compactive  effort  that  is  expected  in  the  field.   The  rela- 
tionship between  the  moisture  and  the  density  is  called  the 
"control  curve".   The  second  stage  involves  the  comparison 
of  field  measurements  of  the  dry  density  with  the  maximum 
dry  density  on  the  control  curve.   Oenerallyi  if  the  percent 
compaction  is  between  957.  and  lOO'At    the  compaction  is  judged 
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to  be  satisfactory.   If  the  percent  compaction  is  con- 
sistently greater  than  1007.,     the  laboratory  test  is  probably 
not  imparting  as  much  compactive  effort  as  the  compactor. 
Conversely!  if  the  percent  compaction  is  consistently  less 
than  95%,  either  the  lab  compaction  is  imparting  more  com- 
pactive  effort  than  the  field  compactor  is  able  to  deliver, 
or  more  passes  of  the  field  compactor  arm   required. 

To  adjust  the  laboratory  compaction  technique  to 
approximate  the  effects  of  the  compactive  effort  in  the 
field,  it  is  helpful  to  quantify  the  value  of  the  compactive 
effort  of  a  desired  compactor  (Selig,  1971).   Selig  proposed 
that  all  compactors  be  represented  by  a  single  drum  toued  by 
a  separate  machine.   In  this  representation,  the  total  com- 
pactive effort  of  the  roller  for  smooth  wheel,  pneumatic 
tire,  and  tamping  compactors  is 

e  «  FLP  <2.1> 

inhere 

e  =  total  compactive  effort  <ft ' lb> 

F  =  towing  force  (lb) 

L  '^    distance  towed  (ft) 

P  =  the  number  of  passes  over  the  distance  L 

This  implies  that  all  of  the  compactive  effort  of  these  com- 
pactors is  provided  by  the  towing  unit  thi-ough  the  drawbar. 
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In  contrast,  the  compactive  effort  of  vibratory  rollers  is 
imparted  by  the  vibratory  energy  of  the  roller  drum.   The 
effect  of  the  touting  unit  on  the  drawbar  is  assumed  to  be 
negligible.   Therefore,  the  total  compactive  effort  of  a 
vibratory  roller  is 


e  =  ZtoT^  <;2.2> 


luhere 


Z   =   net  work  per  cycle  ( f t '  lb /cy c le ) 
Q   =   frequency  of  vibration  Ccyc  les./min> 
T    =   time  to  travel  the  distance  L  ''iT:in> 

Although  equations  2.  1  and  2.2  both  assuTie  that  all  of  the 
mechanical  energy  of  the  roller  is  transformed  into  compac- 
tive effort,  only  a  fraction  of  the  mechanical  energy  is 
actually  transformed.   The  ratio  of  the  compactive  effort  to 
the  mechanical  energy  is  called  the  "efficiency".   The  com- 
pactive effort  per  unit  volume  may  be  expressed  with  quasi- 
analytic  expressions  which  are  specific  to  a  compactor. 
These  expressions  are    given  in  Table  2.  1.   All  these  expres- 
sions contain  a  coefficient  of  rolling  friction,  f.   Like 
compactive  effort,  this  coefficient  has  almost  never  been 
determined  experimentally  (Selig,  1971).   To  overcome  this, 
Selig  estimated  the  compactive  effort  of  several  types  of 
compactors  by  comparing  the  densities  that  they  could 
achieve  with  the  densities  obtained  in  Standard  and  Modified 
Proctor  compaction  tests.   Substituting  these  values  into 
the  expressions  in  Table  2. 1»  he  evaluated  the  values  of  f. 
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Table  2.1  Expressions  for  the  Compactive  Effort 
of  Various  Types  of  Compactors  (after 
Selig,  1971) 


Roller  Type  Compactive  Effort 


Smooth  Wheel  E  =  ^^ 


Pneumat  i  c  E  =  rr 


f  UIP 
E  =  ^  if  d  <  2b 
Et 


_  _  fWiTCd  ^^  2> 
Tampir.Q  t  -   ^^   ^  ^^^ 

o 


fWP 
Et 


375H  P 
_    v_ 

...'ibratory  t  =  ~^^ 


fWP 
Bt 


375H 

V 

inhere  f  '^  ~rjr — 
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Table  2.  1   (continued) 


2 

A   =  contact  arec  of  tamping  foot  <ft  > 

B   =  roller  width  (ft) 

b   =  tire  width  (ft) 

c   =  foot  aTea    correction  factor  >  1.0 

D   =  diaiTieter  of  roller  drum  <ft>  ^ 

E   =  compaction  effort  per  unit  volume  >Cft  '  lb/ft  > 

f   =  coefficient  of  compaction 

H   =  horsepouier  of  vibrator  engine 

h^  =  nb  for  d  >  2b 

=  b  +  (n-l)d    for    d  <  2b 

k   =  overlap  correction  factor  <  1.0 

1   =  tamping  foot  length  <ft> 

N'   =  number  of  tamping  feet 

n   =  number  of  tires 

P   =  number  of  passes 

S   =  foriuard  speed  (miles/hour) 

t   =  compacted  lift  thickness  (ft) 

U   =  total  weight  of  compactor  (lb) 
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Thereforei  Selig's  f  coefficient  is  not  actually  a  coeffi- 
cient of  rolling  resistance^  but  a  general  purpose  correc- 
tion factor  called  the  coefficient  of  conpaction<  which 
reflects  the  number  of  compactor  passes<  the  soil  lift 
thickness<  and  the  soil  type.   The  range  of  values  of  f  and 
the  recommended  average  design  values  are  given  in  Table 
2.  2. 

The  relations  for  compactive  effort  in  Table  2.  1  assume 
that  the  compactive  effort  delivered  to  the  soil  is  linearly 
proportional  to  the  number  of  passes  of  the  conpactor.   In 
practice>  this  is  not  the  case.   As  the  number  of  passes  of 
the  compactor  increases*  the  densif icatior;  per  pass  dimin- 
ishes.  This  implies  that  as  the  number  of  passes  increases< 
a  lesser  amount  of  the  work  done  by  the  compactor  is 
transformed  into  compactive  effort.   Selig  (1971)  showed 
that  this  phenomenon  could  be  incorporated  in  his  model  with 
the  following  expression  for  the  f  parameter: 


f .  =  kt/P  <2.3> 

1 


where 


th 
f .  =*  coefficient  of  compaction  after  the  i 

pass 

k  B  compaction  constant 

t  s  compacted  layers  thickness 

P  =  number  of  passes 

It  follows  that  the  average  value  of  f  over  P  passes  is 
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Tablt  2.2   Valuas  of  f  and  k  for  us*  in  Tablv  2.1 
and  Equation.  2.6  (after  Selig*  1971) 


Roller  Type 

Soil  Type 

f 

0.  35 

k 

Sheepsfoot 

2 

0. 

20-0.  50 

1.9 

3 

0. 

05-0.  15 

0.  10 

1.  6 

Pneumatic 

1 

0. 

05-0.  25 

0.  15 

1.  1 

2 

0. 

lO-O.  30 

0.25 

1.  3 

3 

0. 

05-0.  25 

0.  15 

1.  1 

4 

0. 

05-0.30 

0.  15 

1.  4 

Smooth  Wheel 

1 

0. 

20-0.  50 

0.  35 

1.  3 

2 

0. 

15-0.25 

0.  15 

1.  1 

3 

0. 

20-0.  50 

0.25 

1.  1 

4 

0. 

10-0.40 

0.  30 

1.  3 

Vibratory 

1 

0. 

20-0.  40 

0.  30 

2.  7 

2 

0. 

15-0.  40 

0.  25 

2.  2 

3 

0. 

3O-0.  60 

0.  40 

2.  5 

4 

0. 

50-1.  00 

0.  SO 

2.  7 

Segmented  Pad 

2 

0. 

10-0.30 

0.  20 

0.  9 

3 

0. 

10-0.  25 

0.  15 

1.  1 

Soil  Type       Description 


1  sand 

2  silty  clay 

3  gravel-sand— clay 

4  crushed  stone 
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4> 


u  h  £'  r  e 


f  =  k ■ t  M/P 


M  =   E  - 


Substituting  f  into  the  basic  expression  for  compactive 
effort  (Table  2.1)  yields: 

M 


E  =  k  W  ^ 


'■.  2  . 6  ':• 


Suggested  values  of  k  for  use  in  equation  2.6  are  given  in 
Table  2.2    Unlike  the  equations  in  Table  2  1,  equation  26 
does  not  reflect  the  effect  of  variables  specific  to  indivi- 
dual compactors  such  as  tire  spacing-  tamping  foot  length, 
vibrator  horsepower  and  operation  speed.   This  simplifica- 
tion uas  necessary  due  to  insufficient  data.  Even  so,  the 
compactive  effort  should  be  estimated  with  equation  2.6 
because-  unlike  the  equations  in  Table  2.2/  it  simulates  the 
dimunition  in  compaction  per  pass.   The  t\i}Q    formulations  of 
the  compactive  effort  equations  are  compared  in  the  follow- 
ing two  examples  for  the  compactors  which  were  used  in  the 
Purdue  embankment  performance  study  (Terdich,  19B1). 

Examp  le  2.1 


It  is  desired  to  use  a  Caterpiller  Model  B25  segmented 
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pad  tamping  roller  to  compact  one  foot  lifts  of  St.  Croix 
clay.   The  specifications  for  this  roller  are  provided  in 
Table  2.3.   Determine  the  variation  of  c  o.-:ipac  t  i  ve  effort 
delivered  to  the  soil  with  the  number  of  passes  of  the 
roller. 

The  total  ujeight  of  the  roller  is  approximately  60/000 
pounds.   The  compaction  is  performed  uith  four  drums  each 
with  a  44.5  inch  width.   Using  the  expression  for  compactive 
effort  of  tamping  rollers  given  in  Table  2. 1/  the  average 
compactive  effort  is: 


^  <44.5/12>0>     Ifc-.lt.- 


The  coefficient  f  may  be  assumed  to  be  approximately  0.2 
(using  the  value  for  silty  clay  in  Table  2.2).   It  follows 
that: 

E  =  3236  P    Cft  ■  lb/ft'*::' 

Alternately/  the  compactive  effort  may  be  evaluated 
using  equation  2.6.   The  k  coefficient  rriay  be  approximated 
by  the  value  of  k  for  a  segmented  pad  conpactor  operating  on 
silty  clay,  i.e.,     k  =  0.9.   Substituting  into  equation  2.6 
yields 


p  _  <0.9><60.000><h>  _       , 
^    ~  <44.5./12>      ~  iH...^ 


The  results  of  this  calculation  are    given  in  Table  2.4.   The 
results  of  the  two  expressions  for  compactive  effort  are 
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Table  2.3   Specifications  of  Compactors  Used  in  Purdues 

Embankment  Performance  Project  (after  Terdich, 
1981) 


Caterpillar  Model  825 


Dimensions 

Length,  with  dozer  23  ft,  4  in. 

Width,  w/o  clearers  11  ft,  11  in. 

uj/o  dozer  12  ft,  6    in. 

u)/o  dozer  13  ft,  7  1/2 

Wheelbase  140  in. 

Weight  (shipping) 

w/o  dozer  59,000  lb. 

with  dozer  63,000  lb. 

Number  of  Drums  4 

Number  of  Pads/Drums  65 

Each  drum  width  44  1/2  in. 

Max.  ballast  244  U.  S.  Gal 
Bulldozer  Dimensions 

Length  14  ft. 

Height  40  1/2  in. 

Maximum  Speeds  <?ear      MPH 

1  3.  1 

2  7.  0 

3  17.  0 


Table  2.3   (continued) 


RauQo  Rascal  Model  42QC 
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Dimensions 

Lengthf  uith  blade 

Width 

Wheelbase 


Weight 

Vibration  Drive 

Frequency 

Dynamic  Force 

No.  of  Drums 
No.  of  Pads/Drums 
Each  drum  uiidth 
Maximum  Speeds 


18  ft,  9  in. 
9  ft 
9  ft 


25,  160  lb. 
Hydraulic,  Direct  Drive 
1100  to  ISOO  rpm 
32,000  lb 


1 

140 

84    in. 

Gear 

MPH 

1 

4 

2 

6 

3 

8 

23 


Table    2.4      Calculations    -    Examples    2.1    and    2.2 


Compac  t  i  ve   Effort 
■•'ft  •  lb./ft    > 


Passes 

M 

E  =  14562  M 

E  =  7907  M 

2 

1.  50 

21843 

11861 

4 

2.  08 

30289 

16447 

2 

1.  50 

21843 

11861 

4 

2.  08 

30289 

16447 

6 

2.  45 

35676 

19373 

8 

2.  72 

39608 

21508 

10 

2.  93 

42666 

23169 

12 

3.  10 

45142 

24513 

14 

3.  25 

47326 

25699 

16 

3.  38 

49219 

26727 

2A 


cofflparvd  in  Figure  2.  4.   The  equation  for  compactive  effort 
of  Table  2.2  underestimates  the  compactive  effort  up  to  15 
passes  because  it  does  not  simulate  the  dimunition  in  com- 
paction per  pass.   Note*  that  the  values  of  compactive 
effort  for  this  roller  are  generally  intermediate  to  the 
total  values  of  the  Standard  and  Modified  Proctor  compaction 
tests. 

Example  2.  2 

Repeat  Example  2.  1  for  a  RayGo  Rascal  Model  420C  vibra- 
tory roller.   Specifications  are  provided  in  Table  2.3.   The 
width  of  the  drums  is  84  inches.   The  total  weight  of  the 
compactor  is  25. 160  pounds.   Assume  that  the  roller  proceeds 
at  1.  5  miles  per  hour.   The  f  coefficient  is  estimated  to 
have  the  value  given  in  Table  2.  2  for  a  vibratory  roller  on 
silty  clay>  i.e..  f  =  0.25.   The  horsepower  of  a  vibrator  is 
<Table  2.  1): 


H   =  0.0027  fWS  <2.7> 

V 

Therefore,  the  horsepower  of  this  particular  vibrator  is: 
H   =  <.0027><. 25X21, 600><1.5>  =  21.87  hp 

V 

The  expression  for  the  compactive  effort  of  a  vibratory 
roller  from  Table  2.  1  is: 
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<375><21  .87XP 
<1  .5X84/1 


:,: :'  :  =  78i  p 


Al  ternately  1  the  compactive  effort  iTay  be  obtained  from 
equation  2.6.   Note  that  equation  2.6  employs  the  static 
u/eight  instead  of  the  dynamic  force.   This  discrepancy  is 
accounted  for  in  the  k  coefficient  (k  =  2.2,  see  Table  2.2), 
Substituting  into  equation  2.6  yields; 

=  <2.2><25..160XM:.  ^ 
C84/12> 

The  results  of  this  calculation  are    presented  in  Table  2.4 
and  Figure  2.  5.   Although  equation  2.  6  does  not  include 
speed  of  the  roller*  it  does  account  for  the  reduction  in 
compaction  per  pass  as  compaction  proceeds.   Therefore,  as 
a'as  the  case  uith  the  tamping  roller,  the  compactive  effort 
of  the  vibratory  roller  is  better  predicted  by  equation  2.6 
than  the  expression  in  Table  2. 1.   This  particular  compactor 
will  yield  the  value  of  total  Standard  conpactive  effort 
after  approximately  tiuo  passes  (See  Figure  2  5).   A  com- 
parison of  Figures  2.4  and  2.5  indicates  that  the  Cater- 
piller  tamper  is  expected  to  be  more  effective  than  the  Ras- 
cal Vibratory  roller  for  compacting  clay. 

Assuming  a  selection  of  compactor  and  its  use,  based 
upon  experience,  a  laboratory  test  may  be  selected,  assuming 
that  the  total  energy  in  the  field  and  laboratory  are    equal. 
The  equation  for  calculating  the  compactive  effort  in  a  com- 
paction test  is: 
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E    =    H    W    B    L/y  <2.8> 

(uhere 

E  =  compactive  effort 

H  =  height  of  fall  of  the  harnmer 

W  =  weight  of  the  hammer 

B  =  number  of  blotus  per  soil  layer 

L  =  number  of  soil  layers 

V  =  volume  of  the  mold 

If  equipment  is  already  available  to  perform  an  impact 
compaction  testi  the  volume  of  the  mold,  the  lueight  of  the 
hammer,  and  the  height  of  fall  are  predetermined.    If  the 
number  of  layers  is  set  as  in  the  Standard  and  Modified  com- 
paction tests,  the  only  unspecified  variable  is  the  number 
cif  blou/s  per  layer.   This  may  be  adjusted  as  desired  to 
approximate  the  rollers  compactive  effort.   Rearranging 
equation  2  8.  the  required  number  of  blous  is: 


'..'  E 
B  =  f^-^^j-j-  <2.Sa> 


Examp  le  2.  3 

It  is  desired  to  prepare  a  soil  sample  to  the  level  of 
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compaction  expected  from  the  Rascal  compactor  after  four 
passes  (Example  2.2).   The  sample  will  be  compacted  with  the 
equipment  used  for  performing  standard  compaction  tests: 

V  =  1/30  ft'^ 

L  «=  3  layers 

W  =  5.  5  lb. 

H  =  1.  0  ft. 

From  Table  2.4.  the  energy  corresponding  to  four  passes  is 
16447  ftlb/ft^.   Therefore  the  number  of  blows/layer  should 
be 

^  <1/30><16447>  ^     blows/layer 
<1><5.5><3>         oiows^i«y«r 

This  represents  approximately  25%  more  compactive  effort 
than  standard  compaction. 

The  Hubrid  Approach 

Figure  2.2  shows  that  the  OMC  varies  along  the  line  of 
optimums,  and  as  the  compactive  effort  increases,  the  OMC 
decreases.   When  a  compactor  imparts  a  compactive  effort 
associated  with  an  OMC  equal  to  the  water  content  in  the 
field,  the  compaction  process  is  efficient.   Therefore,  the 
specified  compactive  effort  should  correspond  to  an  OMC 
equal  to  the  expected  compaction  water  content.   This  combi- 
nation of  specification  of  results  and  specification  of  pro- 
cedure involving  a  field  compactive  effort  is  called  the 
hybrid  approach. 
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Example  5.  4 

The  insitu  water  content  of  St.  Croix  clay  that  is  to 
be  used  to  construct  an  embankment  is  20%.   The  variation  of 
the  compactive  effort  vs.  the  OMC  for  laboratory  compaction 
is  Bhoun  in  Figure  2.  h.        The  compaction  ii>ater  content  is 

equal  to  the  OMC  for  a  compactive  effort  of  approximately 

3 
25,000  ft  lb/ft  .   Assuming  that  the  efficiency  of  the  com- 
paction in  the  field  is  approximately  equivalent  to  labora- 
tory compaction,  three  passes  of  Caterpillar  Model  825 
should  be  sufficient  to  impart  this  compactive  effort  to  the 
soil  (Figure  2.4).   Once  the  choice  of  compactor  has  been 
made,  a  density  growth  curve  should  be  developed  on  a  test 
section  to  check  the  prediction. 
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FIGURE  2.6  VARIATION  OF  OPTIMUM  MOISTURE  CONTENT 

VS.  COMPACTIVE  EFFORT  FOR  ST.  CROIX 
CLAY  (AFTER  D  I  BERNARDO.   I  979) 
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III  -  SLOPE  STABILITY  CONSIDERATIONS  FOR  EMBANKMENT  DESIGN 

The  Concept  of  the  Factor  of  Safetu 

The  first  criterion  to  be  satisfied  in  the  design  of  a 
compacted  embankment  is  the  stability  of  the  embankment  side 
slopes.   Typically*  the  relative  stability  of  a  slope  is 
assessed  uith  a  factor  of  safety  obtained   by  limit  equili- 
brium analysis.   The  factor  of  safety  is  defined  as  the 
ratio  of  available  strength  to  applied  shear  stress  along  a 
surface  of  unit  thickness  beneath  the  free  surface  of  the 
slope.   Each  slope  has  a  family  of  such  surfaces.  The  sur- 
face u>ith  the  minimum  factor  of  safety  is  referred  to  as  the 
critical  surface.   If  the  factor  of  safety  on  the  critical 
surface  is  greater  than  one*  the  slope  is  considered  stable. 
Conversely,  if  the  factor  of  safety  on  the  critical  surface 
is  less  than  one,  the  slope  is  considered  unstable.   Since 
the  value  of  the  factor  of  safety  that  controls  a  design  is 
the  minimum  value,  the  minimum  factor  of  safety  will  be 
referred  to  as  the  factor  of  safety  throughout  this  chapter, 
unless  specified  otherwise. 

Factors  that  complicate  the  relationship  between  the 
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critical  surface  and  the  expected  failure  surface  include: 

(1)  The  dependence  of  the  position  of  the  critical 
surface  on  the  factor  of  safety. 

(2)  The  deviations  of  the  soil  shear  strength  behavior 
from  the  mathematical  models  used  to  quantify  it. 

(3)  Errors  inherent  in  the  way  slope  stability  analysis 
methods  calculate  the  normal  stresses  along  the 
trial  surface. 

(4)  Additional  resistance  due  to  end  effects  in  actual 
three-dimensional  surfaces. 

In  structural  mechanics,  the  factor  of  safety  is  com- 
puted by  comparing  the  resistance  of  a  structure  to  a     load 
that  IS  applied  to  the  critical  surface.    For  example/   if  an 
axial  load  is  applied  to  the  prismatic  e lasto-p last i c  bar 
shooin  in  Figure  3.  1,     the  critical  surface  will  be  the 
cross-section  with  the  smallest  area.    Therefore,  the  factor 
of  safety  will  be 


FS  =  a  ■  A   /P  i;^.  1 

y   cr 


u.'  h  e  r  e 


cr   =  yield  strength  of  the  bar 
rt_   =  area  of  the  smallest  cross-section  of  the  bar 
P  =  the  applied  load 


Tne  factor  of  safety  in  this  case  is  a  measure  of  the  prox- 
imity to  failure  of  the  cross-section  on  which  the  bar  is 
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FIGURE     3. 1  CRITICAL  SURFACE  OF  AN  AXIALLY 

LOADED  TAPERED    PRISMATIC  BAR 
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expected  to  fail  if  the  load  is  increased  or  the  yield 
strength  lowered.   This  is  not  the  case  for  the  factor  of 
safety  used  in  slope  stability.   Consider  the  embankment 
shown  in  Figure  3.2a,  whose  geometry  is  defined  by  the  slope 
hieight.,  H..  and  the  slope  angle ^  B.       The  shear  strength  is 
assumed  to  be  a  constant  everywhere  in  and  under  the  embank- 
nent.   If  the  value  of  the  soil  density  increases,  the  fac- 
tor of  safety  will  decrease  and  the  critical  surface  will 
move  progressively  deeper  under  the  embankment  provided  that 
the  soil  strength  does  not  increase  because  of  the  density 
increase  (Figure  3.2b).   Each  value  of  the  factor  of  safety 
corresponds  to  a  different  critical  surface.    This  implies 
that  the  slope  is  expected  to  fail  along  a  different  surface 
than  the  critical  surface. 


Typically,  a  shear  strength  envelope  is  obtained  by 
using  peak  values  of  the  deviatoric  stress  from  a  triaxial 
test  run  to  simulate  insitu  conditions.   However,  the  soil 
in  the  slope  may  only  be  able  to  sustain  a  reduced  devia- 
toric stress  because  of  strain  softening.    Since  the  states 
of  stress  and  strain  vary  greatly  from  position  to  position 
within  an  embankment,  it  is  unlikely  that  the  maximum  values 
of  the  strength  envelope  can  be  developed  simultaneously 
along  any  trial  surface  on  which  the  factor  of  safety  is  to 

hp     o\;5liia^oH  TharoPrtTp.       iiciriu      a      cfTpng*h      pn.'o](ipp     hAc;o({      nn 

peak  values  of  the  deviatoric  stress  will  usually  result  in 
an  ovorsstimate  of  the  factor  of  Baf»ty.    Although  th» 
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FIGURE     3.2a)  SIMPLE  SLOPE 


\>  r^ 


FIGURE    3.2b)      CRITICAL  SURFACES  ON 
A  SIMPLE  SLOPE 
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unconiol idated-undrained  <UU)  shear  strength  envelop*  of 

unsaturated  compacted  clay  is  nonlinear  if  the  water  content 

is  less  than  the  OMC  (Weitzel.  1979).  it  may  be  represented 

by  a  straight-line  Mohi — Coulomb  envelope  (Figure  3.3).   This 

simplification  underestimates  the  strength  along  portions  of 

tfie  surface  uihere  the  normal  stress  <a)  is  iri  the 

CT   <  CT  <  <T.  range  <Zone  II  >  and  overestimate  the  strength 
a        b 

outside  these  limits  (Zone  I).   The  net  effect  can  cause  the 
factor  of  safety  to  be  overestimated  or  underestimated 
depending  on  the  percentage  of  the  trial  surface  that  is  in 
Zone  I  or  Zone  II. 

Generally/  the  parameters  that  define  the  density  and 
strength  of  soil  in  and  under  a  slope  will  vary  with  posi- 
tion.  In  order  to  analyze  such  cases,  it  is  current  prac- 
tice in  slope  stability  analysis  to  resort  to  a  method  of 
slices.   Many  such  methods  exist,  each  with  its  own  set  of 
assumptions  to  overcome  the  indeterminacy  of  the  problem. 
These  methods  assume  that  the  local  factor  of  safety,  i.  e.  . 
the  factor  of  safety  on  a  particular  slice,  is  equal  to  the 
global  factor  of  safety  of  the  soil  mass  above  the  trial 
surface.   This  is  only  exactly  true  at  lifTiit  equilibrium. 
Consider  the  unsaturated  slice  in  Figure  3.4.   Assuming  that 
the  slices  are  of  a  negligible  width,  the  forces  on  the 
sides  of  the  slice  are  equal  and  opposite.  The  summation  of 
forces  in  the  direction  normal  to  the  bottom  of  the  slice 
gives: 
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Zone  I,  in  Mohr  Coulumb  envelope  overestimates    shear  strength 
Zone  n    Mohr  Coulumb  envelope  underestimates  shear  strength 


FIGURE  3.3  COMPARISON  OF  AN  ACTUAL  ENVELOPE 

AND  A  MOHR-COULUMB  REPRESENTATION 
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X-t-AX 


E  +  AE 


ASr=(ANtan</)a  +  ACa)/FS 


FIGURE  3.4 


FORCES  ON  A  SLICE  (AFTER  SIEGEL.   I  975  ) 


itO 


£,U   -    CAW  -  AS   siTia3coso  <3.2) 

where 

AN  =  the  force  normal  to  the  bottom  of  the  slice 
AW  =  the  uieight  of  the  slice 
AS   =  the  resisting  force  along  the  bottooi  of 
the  slice 
a  =  the  angle  that  the  bottom  of  the  slice  makes 
luith  the  horizontal  axis. 

If  a  Mohr-Coulomb  envelope  is  used<  then 


AS   =  CAN  tan-f   +  AC  D/FS  <3.3) 

r  a  9 


where 


tan*   =  slope  of  the  Mohi — Coulomb  envelope 
AC   =  the  product  of  the  arc  length  of  the 

bottom  of  the  slice  and  the  Mohr-Coulomb 
cohesion  intercept 
FS  =  local  factor  of  safetg  of  the  slice. 


Substituting  equation  3.3  into  equation  3.2  and  rearranging 
yields: 

AW'cosa  -  AC  ' s i na ' c osa/FS 

AN  =  — r— rr — T^ — ■ 7^ <3.4> 

1    +    tan*    sina    cosa/FS 

where   C      =   c      dx/cosa   and   c       is   the   Mohi — Coulomb   cohesion 
a  a  a 


intercept. 


1^1 

Sine*  the  v«lue  of  &N  in  equation  3.4  is  dependent  on  a 
value  of  the  local  factor  of  safety  that  is  exact  only  at 
limit  equilibrium^  slope  stability  methods  inherently  calcu- 
late incorrect  values  of  the  normal  stress  on  the  surface  on 
which  the  factor  of  safety  is  computed. 

The  extent  of  the  critical  surface  has  been  assuned  to 
be  infinite  in  the  direction  perpendicular  to  the  slope 
cross-section  (Figure  3.  5).  When  the  lateral  extent  of  the 
slope  is  restricted,  the  critical  surface  mill  arc  upwards 
near  the  lateral  boundaries  of  the  slope  (Figure  3.  6).   An 
approach  for  accounting  for  the  three-dimensional  shape  of 
the  limit  equilibrium  surface  on  the  factor  of  safety  was 
developed  by  Chen  (1761).   His  approach  consists  of  the  fol- 
lowing steps: 

1)  Locate  the  critical  circular  surface  for  a  slope 
using  a  slope  stability  method  such  as  the  Simpli- 
fied Janbu   factor  of  safety. 

2)  Calculate  the  factor  of  safety  on  the  Simplified 
Janbu  critical  surface  using  the  Spencer  method. 
This  is  taken  to  be  the  exact  value  of  the  two- 
dimensional  factor  of  safety. 

3)  Assume  that  the  critical  three-dimensional  surface 
is  an  ellipsoid  attached  to  a  cylinder  (Figure  3.6) 
with  a  cross-section  defined  by  the  two-dimensional 
surface. 


1.   The  reader  should  employ  the  Simplified  Janbu  option 
for  circles  with  caution  (Boutrup<  1977). 
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4)  Using  the  LEMIX  program  (Chen<  1931)/  calculate  the 
factor  of  safety  on  the  three-dimensional  surface. 
This  approach  luill  be  illustrated  with  the  folloming  exam- 
ple: 

Examp  le  3.  1 

It  is  desired  to  assess  the  three-dimensional  factor  of 
safety  of  the  slope  shown  in  Figure  3.7.   Assume  that  the 
critical  three-dimensional  surface  is  an  ellipsoid  attached 
to  a  cylinder  (Figure  3.6).   The  cylinder  is  defined  by  the 
radius  of  the  critical  surface  for  the  Simplified  Janbu  fac- 
tor of  safety  which  is  shown  in  Figure  3.8.   The  Simplified 
Janbu  factor  of  safety  on  this  surface  is  1.28.   The  Spencer 
factor  of  safety  on  the  same  surface  is  1.34.   The  width  of 
the  ellipsoids..  l_j  is  defiried  by  -a  sp-^cified  r.3dius  which 
passes  through  the  endpoint  of  the  cylinder    Using  the  pro- 
gram LEMIX  it  is  possible  to  calculate  the  three-dimensional 
fac+^or  of  safety  for  various  values  of  the  l./H  ratio  where 
H  =  height  of  the  slope 
1_  =  ha  If -width  of  the  cyl  iridrical  portion  of 

the  sliding  mass,  i.e.,  the  portion  which  is 
identical  to  the  two-dimensional  rotational 
surface  (see  Figure  3.9). 

The  results  are  presented  in  Table  3.  1  and  Figure  3.  10.   As 

the  1  /H  ratio  increases.,  the  three-dimen-ional  factor  of 

c 

safety  approaches  the  conventional  two-dimensional  value. 
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FIGURE  3.8 


CRITICAL  CIRCLE  FOUND  WITH  SIMPLIFIED  JANBU 
FACTOR  OF  SAFETY 
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Table  3.  1   Three-Dimensi onal  Factor  of  Safety  vs. 
1  _  /H  rat  i  ci  -  Examp  1  >=■  3.1 


1  /H 
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However,  at  values  of  1  /H  less  than  approximately  3  0,    the 
three-dimensional  factor  of  safety  is  considerablg  higher 
than  the  two-dimensional  value.   This  fact  can  be  taken 
advantage  of  in  embankments  of  small  width  because  the  most 
critical  surfaces  (i.e..  those  that  have  a  long  central 
cylinder)  cannot  develop,  due  to  geometrical  limitations. 

Comments  on  the  STABL  Program 

5TABL  is  a  general  purpose  slope  stability  computer 
program  that  was  developed  for  the  Indiana  Department  of 
Highways.   It  is  designed  to  generate  a  number  of  trial  sur- 
faces that  is  specified  by  the  user  and  to  calculate  the 
factor  of  safety  on  these  surfaces.   The  user  specifies  the 
zone  in  which  he  desires  the  surfaces  to  be  generated.   This 
distinguishes  STABL  from  other  slope  stability  programs  with 
searches  specified  by  central  coordinates  and  radii  of  cir- 

(  1 1  1  A  r   art  • 

The  STABL  program  does  not  actually  find  the  critical 
surface  since  there  is  an  infinite  number  of  possible  random 
shaped  surfaces.   Instead*  STABL  outputs  the  ten  most  criti- 
cal surfaces  that  mre    found  and  their  respective  factors  of 
safety.   These  surfaces  are  plotted  automatically.   An 
advantage  that  results  from  this  methodology  is  that  the 
factor  of  safety  on  various  paths  through  a  slope  can  be 
compared.   For  example,  consider  the  embankment  in  Figure 
3.11.   In  a  complete  stability  analysis,  the  critical 
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surface  may  be  found  to  pass  through  till  and  clay  layers 
that  underly  the  embankment.   The  corresponding  factor  of 
safety  is  an  overall  factor  of  safety.   Hou;ever<  if  the 
search  is  restricted  to  surfaces  that  pass  exclusively 
through  the  embankment  itselfi  the  factor  of  safety  on  the 
most  critical  surface  u/ithin  the  embankment  is  an  "intrin- 
sic" factor  of  safety  that  is  unique  to  the  embankment.   The 
intrinsic  factor  of  safety  is  frequently  an  upper  bound  to 
the  actual  factor  of  safety  because  surfaces  with  lower 
values  of  the  factor  of  safety  may  be  found  through  layers 
that  underly  the  embankment. 

STABL  has  several  different  surface  generation  options. 
A  short  description  of  these  options  is  included  in  the  fol- 
louing  paragraphs. 

The  SURFAC  option  is  a  command  that  may  be  used  to 
determine  the  factor  of  safety  on  a  surface  of  general  shape 
specified  by  the  user.   The  factor  of  safety  on  the  surface 
is  ralnilatprt  by  th»  Rimplifiod  v'anbii  mothnd     Tho  SURTAC 
option  is  used  to  check  the  factor  of  safety  calculated  by 
SIAQL  against  documented  solutions. 

The  SURBIS  option  is  identical  to  the  SURFAC  option 
uith  the  exception  that  it  computes  the  factor  of  safety 
with  the  Simplified  Bishop  method.  Care  should  be  taken  that 
a  circular  surface  is  input  because  the  Simplified  Bishop 
method  calculates  an  incorrect  value  of  the  factor  of  safety 
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if  th«  coo1*dinatc«  of  th»  limit  vquilibriufli  surfact  art  not 
circular  (Bishop.  1995). 

The  CIRCLE  option  randomly  generates  circular  surfaces 
and  evaluates  the  corresponding  factors  of  safety  uith  the 
Simplified  Janbu  factor  of  safety.  Searching  with  circular 
surfaces  generally  yields  a  critical  surface  whose  factor  of 
safety  is  nearly  as  low  as  the  factor  of  safety  found  on  the 
critical  noncircular  surface. 

The  CIRCL2  option  is  identical  to  the  CIRCLE  option 
except  that  the  factor  of  safety  is  calculated  by  the  Sim- 
plified Bishop  method.   This  option  is  generally  preferred 
to  the  CIRCLE  option  because  there  is  more  experience  with 
the  Simplified  Bishop  method  than  with  the  Simplified  Janbu 
mCfthod   Also>  the  Simplified  Biehop  method  does  not  h«v(B  thp 
c  onv«>r  y  enc  e  problemtt  som*tlm«b  encountered  with  the  Simpli- 
fied Janbu  method  for  slopes  that  have  high  cohesion  inter- 
cepts and  low  friction  angles. 

In  nonhomogeneous  slopesi  the  critical  surface  may  be 
noncircular.   The  same  situation  may  arise  if  a  homogeneous 
slope  is  subjected  to  a  surcharge  load  or  pseudo-static 
earthquake  loading.   STABL  supplies  the  RANDOM  option  for 
these  cases.   The  RANDOM  option  pseudo-randomly  generates 
noncircular  surfaces  and  evaluates  the  factor  of  safety  on 
them  with  the  Simplified  Janbu  method.   An  example  of  a  ran- 
domly generated  noncircular  surface  is  shown  in  Figure  3. 12. 
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The  BLOCK  option  specifies  a  straight  line  surface 
bctuieen  randomlg  chosen  points  in  boxes  of  a  size  specified 
by  the  user  (Figure  3. 13).   The  remaining  portions  of  these 
trial  surfaces  are  generated  randomlg  to  the  left  of  the 
leftmost  box  and  to  the  right  of  the  rightmost  box.   The 
Simplified  Janbu  method  is  used  to  calculate  the  factor  of 
safety  on  these  surfaces.   The  BLOCK  option  is  especially 
effective  when  there  is  a  weak  seam  or  bedding  plane  in  the 
slope. 

Using  a  method  of  slices  for  the  analysis  of  block 
shaped  surfaces  is  consistent  with  the  factor  of  safety  used 
for  surfaces  of  other  shapes  because  it  imposes  the  condi- 

tinn     that:     fhp      Inisl      f«rtnr     nf     sat  0  t\i     and      t:he     gluhal      faiitii- 

of  safety  are  equal  everywhere  along  the  trial  surface. 
This  is  different  than  ordinary  methods  of  calculating  the 
factor  of  safety  of  a  block  shaped  surface.   These  methods 
implicitly  assume  that  the  local  factor  of  safety  on  the 
passive  wedge  and  the  active  wedge  is  unity  although  the 
factor  of  safety  on  the  central  block  must  be  greater  than 
unity  for  stability. 

The  BL0CK2  option  is  identical  to  the  BLOCK  option  with 
the  exception  that  the  portions  of  the  trial  surface  to  the 
left  of  the  leftmost  box  and  to  the  right  of  the  rightmost 
box  are  generated  to  simulate  Rankine  active  and  passive 
zones  respectively.   This  option  was  developed  because  the 
BLOCK  option  yields  a  minimum  factor  of  safety  on  surfaces 
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that  approximate  the  Rankine  state.   Therefore,  use  of  the 
BL0CK2  option  pprmits  a  savirigs  in  computational  pffort 
Mit  ULUCK  option  h<as>  been  letamedi   hotDevei-.   for  use  in 
h&terogeneous  soil  profiles. 

While  providing  advice  on  the  use  of  STABL  in  the  capa- 
city of  STABL  consultant  at  Purduei  the  author  noted  that 
certain  problems  arose  repeatedly.   These  problems  included: 

(1)  surfaces  u/ith  unacceptable  shapes 

(2)  the  direction  of  surface  generation 

(3)  direction  limits  on  benched  slopes 

(4)  variability  of  results  due  to  the  random 
generation  of  surfaces 

(5)  a  lack  of  experience  uiith  the  Simplified 
Janbu  factor  of  safety. 

Tiiese  problems  are  discussed  in  the  following  paragraphs. 

Occasionally/  when  the  RANDOM  option  is  usedi  one  or 
more  of  the  ten  most  critical  surfaces  that  are  generated 
iiiill  be  concave  or  even  have  a  reverse  curvature  as  shoiurt  in 
Figure  3. 14.   Such  surfaces  are  k inemat ical 1 y  impossible. 
Therefore,  the  factors  of  safety  computed  on  these  surfaces 
should  be  disregarded. 

STABL  assumes  that  a  slope  rises  from  left  to  right. 
Therefore,  it  generates  surfaces  that  progress  from  left  to 
right  (Figure  3. 15a).   If  the  user  attempts  to  input  a  slope 
that  rises  from  right  to  left,  STABL  would  continue  to  generate 
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LTL 


RTL 


Termination  range 


LIL  leftmost  initation  limit 

RIL  rightmost  initation  limit 

LTL  leftmost  termination  limit 

RTL  rightmost  termination  limit 


FIGURE   3.15  a       EXAMPLE  OF  SLOPE  INPUT  BACKWARDS 
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trial  surfaces  from  left  to  right>  which  never  reach 
the  termination  range.   In  this  case*  STABL  outputs  error 
message  RC-06  and  halts  execution  of  the  problem.   This 
error  can  be  avoided  fiimplg  by  inputting  the  data  co  that 
the  felope  rlbeft  frbm  left  to  right  (Figure  3.  ISb). 

When  a  user  detlrts  the  initiation  limits  to  straddle  a 
break  in  the  ground  surface  where  the  inclination  of  the 
slope  decreases  from  left  to  right,  the  direction  limits 
must  be  compatible  with  each  of  the  segments  that  lie 
between  the  leftmost  initiation  limit  and  the  rightmost  ini- 
tiation limit.   Otherwise,  a  situation  may  develop  where 
STABL  will  generate  a  trial  surface  that  goes  outside  the 
slope  (see  Figure  3.  16).   If  the  trial  surface  does  not 
cross  the  ground  surface  between  the  termination  limits, 
error  message  RC-10  is  output  and  execution  of  the  problem 
is  halted.   If  the  trial  surface  does  cross  the  ground  sur- 
face between  the  termination  limits,  the  results  are  mean- 
ingless, although  no  error  is  detected  by  the  program.   To 
avoid  this  difficulty,  the  user  should  define  left  and  right 
initiation  points  for  each  segment  of  the  ground  surface  in 
the  initiation  range.   Each  set  of  initiation  limits  should 
be  executed  with  appropriate  direction  limits  on  a  separate 
run.   This  is  illustrated  in  Figure  3.17. 

The  inclination  of  the  first  two  segments  of  circular 
and  random  shaped  surfaces  generated  by  STABL  and  the  incli- 
nation of  subsequent  segments  of  random  shaped  surfaces  are 
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FIGURE    3.15  b 


EXAMPLE  OF  CORRECTLY  INPUT  SLOPE 
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pssudo-randomlg  generated  uith  the  aid  of  a  random  number 
generator  called  RANF.   RANF  is  a  computer  supplied  function 
on  CDC  6000  series  computers  that  generates  uniformly  dis- 
tributed random  numbers  between  0  and  1.   However/  since  the 
generation  sequence  starts  with  the  same  built-in  seed  each 
time  that  the  program  is  used<  it  always  yields  identical 
random  number  sequences.   Therefore.  STABL  creates  identical 
surfaces  each  time  it  is  run  on  a  CDC  co-nputer  unless  the 
seed  is  modified.   When  STABL  is  run  on  an  IBM  computer>  the 
user  must  supply  a  random  number  generator.   Some  users  have 
employed  generators  that  yield  a  different  sequence  each 
time  the  program  is  run.   When  this  is  the  case<  STABL  will 
create  different  trial  surfaces  each  ti.T.«?  the  program  is 
runi  and  consequently,  slightly  different  values  of  the  fac- 
tor of  safety. 

Originally,  all  the  STABL  options  for  calculating  the 
factor  of  safety  on  noncircular  surfaces  employed  the  Carter 
method  (Carter,  1971).   This  method  makes  the  following 
assumptions: 

1>  All  interslice  forces  are  equal  in  magnitude  and 
opposite  in  sign.   Therefore,  they  may  be  neglected 
in  the  analysis  of  the  factor  of  safety. 
2)  The  minimum  factor  of  safety  on  any  trial  surface 
is  obtained  when  moments  of  the  driving  and  resist- 
ing forces  are    taken  about  a  point  that  is  at  an 
infinite  height  above  the  slope. 
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3)  All  other  assumptions  are  identical  to  those  made 
bg  the  Simplified  Bishop  method. 
The  expression  for  the  factor  of  safety  given  by  the  Carter 
fltethod  does  not  require  the  calculation  of  interslice 
forces.   Therefore*  the  reasonableness  of  the  line  of  thrust 
need  not  be  checked  as  is  required  by  methods  such  as  the 
extended  Spencer  method  (Spenceri  1773)  and  the  Rigorous 
Bishop  method  (Bishop*  1955).   This  permits  the  Carter  fac- 
tor of  safety  to  be  calculated  with  relatively  little  compu- 
tational effort. 

Boutrup  (1977)  indicated  that  the  Carter  method  is 
identical  to  the  Rigorous  Janbu  method  (Janbu>  1954)  with 
the  simplifying  assumption  that  interslice  forces  may  be 
neglected  from  the  formulation.   This  implies  that  the  Car- 
ter method  (or  the  Simplified  Janbu  method)  does  not  satisfy 
the  requirements  of  statics  for  each  slice*  although  hor- 
izontal force  equilibrium  is  satisfied  for  the  soil  mass 
above  a  trial  surface.   Experience  with  other  incomplete 
equilibrium  formulations  such  as  the  Fellenius  method  of 
slices  indicates  that  incomplete  equilibrium  techniques  give 
values  of  the  factor  of  safety  that  are  lower  than  those 
found  by  complete  equilibrium  techniques.  Since  the  Simpli- 
fied Janbu  method  is  conservative  compared  to  more  exact 
methods*  it  may  be  worthwhile  to  adjust  it  to  correspond  to 
methods  that  yield  results  that  are    closet-  to  complete 
equilibrium  solutions. 
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In  order  to  provide  a  guide  for  performing  this  adjust- 
ment/ values  of  the  factor  of  safety  have  been  computed  for 
simple  homogeneous  slopes  uith  various  values  of  the 
sideslope^  B,    and  the  parameter,  X  =  tant/Cc/VH]  for  both 
the  Friction  Circle  method  and  the  Simplified  Bishop  method. 
The  differences  (in  percent)  of  the  factor  of  safety  of 
these  tuo  methods  relative  to  the  Simplified  Janbu  factor  of 
safety  using  circular  surfaces  uere  also  computed.   The 
results  are  given  in  Table  3.2  and  Figures  3.  18  and  3.  19. 
The  Simplified  Janbu  factor  of  safety  may  be  adjusted  to 
simulate  these  familiar  methods  of  determining  the  factor  of 
safety  uiith  the  following  expression: 


FS  =  FSg^^gj^Cl.O  +  y.   error/lOO:  <:3.5> 

uhere 


FS  =  the  value  of  the  factor  of  safety  that  has 
been  adjusted  to  simulate  another  method. 
FS__^p.  =  the  minimum  factor  of  safety  obtained 
with  the  Simplified  Janbu  irethod 
coded  in  the  STABL  program. 

v  „^^ ^^< desired  met hod >  "  ^^STABL    ,^^ 

/.   error  =  r^ x  100 

STABL 

Values  of  '/.  error  may  be  interpolated  fron  Figures  3.  18  and 
3.  19  for  the  Friction  Circle  and  the  Simplified  Bishop 
methods/  respectively.   This  adjustment  was  developed  for 
simple/  unsaturated/  homogeneous  slopes  with  circular  trial 
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Table  3.2   Factor  of  Safety  for  Values  of  Sideslope  and 
X.   < after  Boutrup,  1977> 


cot  j6 


Met  hod      X        1.5        2.5       3.5 


/.   error 


20 


50 


1  2  1.147  1.443  1.714 

2  1.103      1.351     1.595 

i-rror  3.99       6.81      7.46 


1 

.785 

2 

.409 

3 

.008 

1 

.  703 

2 

.277 

2 

.850 

4 

.82 

5 

.80 

5 

.54 

2 

.368 

3 

.309 

4 

.  226 

2 

.262 

3. 

.  154 

4. 

045 

4 

.69 

4, 

.91 

4. 

.47 

2 

.257 

3, 

.347 

4 

.425 

2, 

.  167 

3 

.242 

4 

.307 

3 

.72 

3 

.24 

2 

.74 

1 

.902 

2 

.948 

3 

.991 

1 

.855 

2 

.893 

3 

.928 

y.   error  2 .  53       1  .  90      1  .  60 


Method  1  -  Friction  Circle 

Method  2  -  STABL2,  Simplified  Janbu 

X=tan-*-/Cc/yH>: 
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surfaces.   However^  it  mag  be  used  approximately  on  slopes 
that  do  not  satisfy  these  conditions. 

It  should  be  emphasized  that  this  adjustment  is  onli/  an 
approximate  procedure  to  correct  an  inco(?.plete  equilibrium 
techniquei  i.e./  a  method  that  does  not  satisfy  equilibrium 
of  the  slices  above  the  trial  limit  equilibrium  surface.   It 
uould  be  more  consistent  to  use  a  method  of  calculating  the 
factor  of  safety  that  satisfies  equilibrium  of  each  of  the 
slices. 

Corrections  to  5TABL 

Originally/  the  STABL  program  calculated  the  Simplified 
Bishop  factor  of  safety  by  multiplying  the  term<?  of  the  fol- 
lowing summation  uhich  is  used  for  calculating  the  Simpli- 
fied Janbu  factor  of  safety 

n    A   +  FS  A 

^   "^  FS  -K  ft   ^  =  ^'  ^^-^'^ 

by   y   =   Rcosa       <see   Boutrup,     1977> 

u'here 

y  =  vertical  distance  to  the  bottom  of  the  slice 

from  the  moment  center 
R  =  radius  of  the  trial  circle 
CT  =  slope  of  the  bottom  of  the  slice 
n  =  number  of  slices 
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A  ,A^^A^  =  terms  reflecting  the  cortdition- 
of  slice  (see  Siegel,  1975). 
FS  =  the  factor  of  safety 

Since  the  value  of  R  is  constant  for  circular  surfaces*  this 
reduces  to: 

n         A   +  FSA 

S   Ccosa  ic   *   £, — ^    =  ^'  <3.7> 

i=l         ""^    "3 

This  procedure  is  very  efficient  for  computer  coding  because 
the  terms  of  the  summation  for  the  Simplified  Bishop  method 
may  be  obtained  by  multiplying  the  terms  of  the  Simplified 
Janbu  summation  by  their  respective  valurs  of  cos  a. 

Unfortunately*  this  formulation  assur.es  that  all  forces 
acting  on  a  slice  act  along  the  base  of  the  slice  <Hou>landi 
1932).   This  is  not  true  if  pseudo-static  earthquake  forces 
or  boundary  surcharge  forces  act  on  a  slice  or  if  the  water 
table  extends  above  a  slice.   Therefore*  STABL  gave 
incorrect  values  of  the  Simplified  Bishop  factor  of  safety 
in  these  circumstances.   To  rectify  this  problem*  the  author 
recoded  the  Simplified  Bishop  factor  of  safety  to  include 
ths  differences  in  the  moment  arms  of  these  forces.   The 
expression  for  the  factor  of  safety  is  (details  of  the 
derivation  are    given  in  Appendix  A): 


n     A. 

f'^l+A-ZFS^ 
FS  = <3.8> 


2A^  -  SA .  +  SA_ 
1  ^    1  ^    1  ^ 
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where 


FS   =    factor    of    safety 

A      =    C    +taTi*'    seca<AW<l-l<    )+iQcos6+6U^+cosje-iU   cosa>    <3.9a> 
1  a  a  V  e  a 

(\      =   tan+'^tana  <3.9b> 

A^   =    CiW<l-k    >    +   iU^   cosjS   +   6Q   cos6>    sina  <;3.9c> 

3  V  .e 

A^    =    (iU^sinjS   +   iQ   sin^Xcosa   -   h/R>  <:3.9d> 

4  jB 

A^   =   k^AQ    Ccosa   -   h      /R>  <3.9e> 

5  h  eq 

These  A-terms  are    not  the  same  as  the  terms  used  for  deter- 
mination of  the  Simplified  Janbu  factor  of  safety,  although 
the  variables  in  equations  3.  9a  through  3. 9e  are  identical 
to  the  variables  used  when  the  program  was  originally 
developed  <Siegel.  1975).   The  factor  of  safety  must  be  cal- 
culated with  an  iterative  procedure  because  equation  3.8  is 
implicit.   In  order  to  make  the  solution  technique  for  the 
Simplified  Janbu  method  analogous  to  that  of  the  Simplified 
Bishop  method/  the  Newton-Raphson  method  which  was  used  for 
the  Simplified  Janbu  factor  of  safety  (Siegel.  1975)  was 
replaced  with  an  implicit  iterative  technique.   The  implicit 
formulation  of  the  factor  of  safety  is  believed  to  be  more 
efficient  in  achieving  a  convergence  on  the  value  of  the 
factor  of  safety  when  incomplete  equilibrium  analyses  are 
performed*  especially  when  the  factor  of  safety  is  much 
greater  or  much  less  than  the  value  that  is  initially 
assumed  (Boutrup.  1977). 
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Convergence  of  the  iteration  scheme  is  achieved  if  the 
assumed  and  back-calculated  values  of  the  factor  of  safety 
on  a  surface  differ  by  less  than  0.005.   The  maximum  number 
of  iterations  is  limited  to  ten.   If  the  solution  does  not 
converge!  the  coordinates  of  the  surface  and  the  last  back- 
calculated  value  of  the  factor  of  safety  are    output.   The 
editing  that  is  required  to  incorporate  these  changes  into 
the  STABL  code  is  listed  in  Appendix  A. 

Strength  Parameters  of  Compacted  Claus 

The  As-Constructed  Condition 

To  replicate  the  condition  existing  at  the  end  of  con- 
struction! a  clay  sample  compacted  to  simulate  field  compac- 
tive  effort  must  be  tested  according  to  UU  procedures.   This 
means  that  the  sample  is  loaded  quickly  so  that  there  is  no 
time  for  the  excess  pore  pressure  induced  by  the  loading  to 
dissipate. 

Ideallyi  the  loading  should  be  performed  to  simulate 
the  stress  path  that  the  soil  undergoes  in  the  field.   Alsoi 
the  loading  should  reflect  the  rotation  in  the  direction  of 
the  stresses  that  occurs  in  the  embankment.   Unfortunately < 
the  initial  stresses  and  stress  changes  in  an  embankment  are 
not  amenable  to  calculation.   Even  if  this  were  possiblei  it 
uould  be  necessary  to  perform  a  prohibitive  number  of  tests 
along  various  stress  paths  in  order  to  model  the  strength  of 
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the  soil  as  a  function  of  position  in  the  embankment.   Con- 
sequently<  the  current  practice  is  to  simulate  the  \J\J    shear 
strength  of  unsaturated'  compacted  soil  under  an  ordinary 
triaxial  loading.   It  has  been  shouin  that  unsaturated'  com- 
pacted clays  tested  in  this  fashion  have  a  strength  line 
that  is  uniquely  defined  by  their  water  content  and  the  com- 
pactive  effort  they  have  been  subjected  to  (Ueitzelf  1979). 
The  shear  strength  data  for  St.  Croix  clay  are  given  in 
Table  3.3  and  in  Figure  3.20.   These  data  may  be  represented 
as  linear  strength  lines  by  performing  linear  regressions  on 
the  data  in  Table  3.  3.  The  Mohr-Coulomb  parameters  may  be 
obtained  -from  the  strength  lines  through  simple  tri- 
gonometric relationships  (Holtz  and  Kovacs>  1981).   The 
results  of  the  regression  are    given  in  Table  3.4  and  Figures 
3.21  and  3.22.  The  correlation  coefficient'  r«  decreases  as 
the  water  content  of  the  soil  increases.   Even  so>  the 
correlations  are   relatively  high  except  when  the  water  con- 
tent of  the  soil  is  well  above  the  OMC.   Using  Figure  3.21 
and  3.22  it  is  possible  to  estimate  the  UU  Mohr-Coulomb 
strength  parameters  of  St.  Croix  clay  at  any  water  content 
for  the  compaction  level  at  which  the  soil  samples  were 
prepared. 

Example  3.  2 

It  is  desired  to  calculate  the  intrinsic  factor  of 
safety  of  the  St.  Croix  clay  embankment  shown  in  Figure 
3. 23.   The  water  content  of  the  soil  is  19%  and  the 
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Table  3.3   UU  Shear  Strength  of  Compacted  St. 
(after  Weitzel,  1979) 


Croix  Clay 


Compac  t  ion 


b)'/. 


(kPa) 


(  kPa) 


^^1^^3V- 


(  kPa) 


Low  Energy*    20.  75 


132 

0 

132 

115 

0 

115 

208 

133 

346 

229 

276 

505 

283 

276 

559 

249 

414 

663 

322 

414 

736 

91 

0 

91 

104 

0 

104 

179 

158 

317 

197 

133 

335 

192 

276 

468 

226 

414 

640 

97 

0 

97 

79 

0 

79 

116 

158 

256 

139 

158 

277 

148 

276 

424 

135 

414 

549 

133 

414 

547 

56 

0 

56 

52 

0 

52 

73 

133 

211 

59 

138 

197 

73 

276 

349 

83 

276 

359 

66 

414 

480 

67 

414 

471 

>  Lou)  Energy  compaction  has  60'/.  of  the  conpactive 
effort  of  standard  compaction 


Table  3.3   (continued) 
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Compac t i  on 


Standard 


u'/. 


19 


20.  1 


21.  6 


25.  0 


'^^1-'='3V''2 

•^3 

^^1*^3^^ 

(kPa) 

(kPa> 

(kPA) 

164 

0 

164 

146 

0 

146 

262 

133 

400 

294 

133 

432 

381 

276 

657 

422 

414 

836 

212 

0 

212 

191 

0 

191 

308 

13H 

446 

320 

276 

596 

352 

276 

628 

372 

414 

786 

200 

0 

200 

158 

0 

158 

188 

133 

326 

237 

276 

513 

250 

276 

526 

243 

414 

657 

240 

414 

654 

74 

0 

74 

73 

0 

73 

115 

153 

253 

87 

133 

225 

89 

276 

365 

107 

276 

383 

150 

414 

564 

99 

414 

513 

Table    3.3       (continued) 
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Compaction 


Modified 


ui'/. 


13.  8 


15.  0 


16.  3 


19.  0 


^^r'='3-V''2 

■^3 

«T^-a3>^/2 

(kPa) 

(kPa) 

(kPa) 

869 

133 

1007 

828 

133 

966 

1007 

276 

1283 

942 

276 

1218 

1044 

414 

1458 

661 

0 

661 

716 

0 

716 

736 

13a 

874 

885 

276 

1161 

977 

276 

1253 

1033 

414 

1447 

915 

414 

1329 

608 

0 

608 

667 

0 

667 

720 

158 

858 

744 

133 

882 

866 

276 

1142 

774 

276 

1050 

747 

414 

1161 

449 

0 

449 

497 

133 

635 

417 

133 

555 

460 

276 

736 

455 

276 

731 

479 

414 

893 
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Table  3.4   UU  Mohr-Coulomb  Parameters  of  Unsaturated  St. 
Croix  Clay 


i 

\- 

c 

Compact  i  on 

OMC 

w-/. 

(  degrees ) 

(kPa) 

r 

Lotu  Energy* 

24.  0 

20. 

75 

17. 

1 

95. 

8 

.  96 

22. 

0 

14. 

0 

S3. 

6 

.  94 

24. 

0 

5. 

9 

83. 

8 

.  81 

27. 

0 

2. 

2 

55. 

9 

.  62 

Standard 

21.  6 

19. 

0 

24. 

1 

108. 

5 

.  99 

20. 

1 

17. 

7 

155. 

1 

.  98 

21. 

6 

B. 

8 

152. 

4 

.  90 

25. 

0 

6. 

1 

67, 

0 

.  78 

Mod  if ied 

16.  3 

13. 

8 

26. 

0 

465. 

4 

.  98 

15. 

0 

26. 

5 

416. 

1 

.  97 

16. 

3 

19. 

0 

461. 

8 

.  88 

19. 

0 

4. 

7 

405. 

8 

.  47 

•iLow  Energy  compaction  has  60'/.  of  the  compactive  effort  of 
Standard  compaction. 
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compactive  effort  is  Standard.   It  can  be  seen  from  Figures 
3.21  and  3.22  that  ♦  and  c  are  24.1°  and  108.5  kPa  <:2266 

psf)  respectively.   The  moist  density  of  St.  Croix  clay  at 

3 
this  uiater  content  and  compactive  effort  is  19.04  kN/m   fl21 

pcf).   Using  STABL*  the  minimum  factor  of  safety  obtained  on 

a  random  shaped  surface  through  the  embankment  is  2.  96. 


The  Long  Term  Condition 

To  replicate  the  behavior  of  an  embankment  long  after 
it  has  been  constructed,  it  is  necessary  to  adjust  strength 
test  procedures.   Long  term  conditions  are  bounded  by  two 
conditions.   The  first  extreme  is  that  long  after  construc- 
tion the  embankment  remains  unsaturated.   The  proper  uay  to 
run  a  test  on  the  soil  for  this  situation  is  the 
consolidated-drained  <CD)  test.   The  soil  is  consolidated  to 
the  expected  state  of  stress  in  the  embankment  and  then 
sheared  at  a  rate  that  is  sufficiently  sloui  to  allou  any 
excess  pore  pressure  developed  by  the  loading  to  dissipate. 
Ordinarilyi  the  requirement  that  the  soil  should  be  consoli- 
dated to  the  expected  state  of  stress  is  relaxed.   Instead, 
the  soil  is  isotrop ical ly  consolidated  to  several  different 
confining  pressures  and  subsequently  sheared.   Experience 
indicates  that  this  has  little  effect  on  the  strength 
envelope  that  is  obtained  for  many  soils.   The  CD  test  is 
not  performed  on  a  routine  basis  because  it  must  be  run  very 
slobJly  to  allow  excess  pore  pressures  to  dissipate. 
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The  second  extreme  condition  is  that  of  the  embankment 
becoming  completely  saturated  sometime  after  construction. 
This  casei  like  the  case  of  the  sudden  drawdown  of  water 
behind  a  dami  requires  the  use  of  a  consol idated-undrained 
(CU>  test.   The  soil  is  consolidated  to  the  expected  state 
of  stress  in  the  field  and  then  sheared  quickly.   As  the 
soil  is  sheared/  the  excess  pore  pressure  due  to  the  stress 
changes  is  measured.   As  was  the  case  with  CD  parameters! 
ths  requirement  of  consolidating  the  sair.ple  to  the  insitu 
state  of  stress  is  usually  relaxed. 

A  major  problem  involved  in  analysis  of  the  second 
extreme  condition  is  the  development  of  excess  pore  pres- 
sure.  The  excess  pore  pressure  that  develops  in  the  labora- 
tory is  caused  by  the  changes  in  the  stresses  that  are 
applied  to  the  sample.   The  excess  pore  pressure  in  the 
field  is  caused  by  the  increase  of  density  due  to  saturation 
and  the  gradual  stress  changes  caused  by  displacements 
within  the  embankment  that  arise  from  changes  from  the  UU  to 
CU  or  CD  conditions.   Therefore*  the  excess  pore  pressures 
measured  in  a  triaxial  test  cannot  be  used  to  predict  excess 
pore  pressures  in  an  embankment. 

An  approximate  method  can  be  used  to  insure  that  the 
excess  pore  pressures  will  be  zero  or  negative.    It  has  been 
shown  that  the  pore  pressure  parameter.,  A-..,  which  relates 
excess  pore  pressure  and  deviator  stress*  varies  with  the 
overconsol idat ion  ratio  (Henkel*  1956).   Typical  results  are 
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shouin  in  Figure  3.24.  The  overconsolidation  ratio  (OCR)  of  a 
compacted  clay  may  be  estimated  uith  the  ratio  of  the  com- 
pactive  prestress  to  the  overburden  pressure.   Therefore, 
the  depth  in  the  embankment  above  which  positive  pore  pres- 
sure can  not  develop  regardless  of  stress  change  (i.e.. 
A^  <  0>  is: 

o 
uihere 

H   =   depth   beneath   embankment    surface   to 

uihich   u   <   0 

u  =  excess  pore  pressure  due  to  the  stress  change 

P   =  c  ompac  t  i  ve  prest  r ess 

y   =  moist  density  of  the  soil 
m  ^ 

OCR   =  the  OCR  at  uihich  u  <  0 
o  — 

Results  in  Figure  3.24  indicate  that  OCR   mill  vary  between 

4  and  5  for  a  natural  saturated  clay.   It  is  reasonable  to 

assume  that  the  same  is  true  for  compacted  clays.   Typical 

results  for  an  embankment  built  of  compacted  St.  Croix  clay 

are  given  in  Table  3.5.   These  results  assume  that 

OCR   =  4.0.   Values  of  P   and  V  are   taken  from  Di Bernardo 
o  s       m 

<1979>. 

Actually,  the  excess  pore  pressure  caused  by  shear  of 
compacted  clays  at  strains  smaller  than  failure  strains  can 
be  larger  than  is  predicted  by  use  of  the  pore  pressure 
parameter  in  Figure  3.24.   Therefore,  OCR   will  not 
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Table  3.5   Depth  to  Zone  of  Positive  Pore  Pressure 
Deve 1 opment 


Compac  t  i  v.-e 

Effort    one  c :.;  ;:■ 


standard 


Modified 


16.  5 


ui  c: '/.  y 

p^  C  k  F'a  > 

H* 

1^  met  ers  ^ 

19.  5 

430.  6 

18.  70 

5.  76 

22.  0 

271.  8 

19.  51 

3.  48 

25.  0 

68.  6 

19.  42 

.  88 

14 

1166.  3 

20.  17 

14.  47 

16.  5 

717.  9 

20.  97 

8.  56 

19.  5 

212.  7 

20.  63 

2.  58 

*H  =  depth  to  the  bottom  of  the  zone  of  positive 
pore  pressure  development 


correspond  exactly  to  a  state  of  zero  excess  pore  pressure 
generation.   In  any  event*  the  depth  above  which  the  excess 
pore  pressure  is  taken  to  be  less  than  zero  mill  sometimes 
be  only  a  small  portion  of  the  embankment  height.   There- 
fore, stress  changes  will  cause  an  increase  in  pore  pres- 
sures.  Unfortunately,  the  stress  changes  and  hence  the  pore 
pressure  changes  can  not  be  predicted  easily.   Conse- 
quently. >  the  only  convenient  approach  is  to  assume  that 
pore  pressure  in  an  embankment  is  equal  to  the  head  of  water 
above  the  location  in  question. 

The  CU  friction  angle  and  cohesion  intercept  of  labora- 
tory compacted  and  saturated  St.  Croix  clay  are   approxi- 
mately 20°  and  15  kPa  regardless  of  compaction  variables  or 
stress  path  <John6on>  1979).   Therefore,  it  is  possible  to 
develop  a  graph  of  the  intrinsic  factor  of  safety  for  every 
possible  geometry.   For  the  unsaturated  case,  the  moist  den- 
sity of  the  embankment  soil  is  used.   For  the  saturated 
case,  the  saturated  density  of  the  embankment  soil  is  used 
and  a  water  table  is  assumed  to  run  along  the  free  surface 
of  the  embankment.   In  these  examples  the  effective  stress 
parameters  are    assumed  to  be  the  same  for  both  saturated  and 
unsaturated  compact-ed  clay.   The  results  are    given  in  Table 
3.6  and  Figures  3.25  and  3.26.   These  results  were  developed 
using  the   Simplified  Janbu  factor  of  safety  on  circular  sur- 
faces that  were  restricted  to  remain  above  the  elevation  of 
the  toe  and  to  exit  the  slope  less  than  10  meters  from  the 
crest  of  the  slope. 
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Table  3.6   Intrinsic  Simplified  Janbu  Factors  of  Safety 
of  Compacted  St.  Croix  Clay  E-nbankments  Using 
Laboratory  Compacted  CU  Shear  Strength 


iJTisaturated  Case 


H  je  FS 


30  2.12 

35  1 . 86 

40  1.71 

45  1.61 

50  1 . 52 

55  1.47 

60  1 . 37 

65  1 . 30 

70  1 . 24 


1  Om  20*^'  2.12 

25  1 . 70 

30  1 . 44 

35  1 . 26 

40  1.1 6 

45  1 . OS 

50  1 . 02 

55  0 . 95 


JOm  17.5  1 . 83 

20 . 0  1 . 59 

22 . 5  1.43 

25 . 0  1 . 29 

30 . 0  1 . 08 

32.5  1.01 


30m  15.0  1.90 

17.5  1.62 

20 . 0  1 . 42 

25.0  1.14 

30 .0  0 . 96 


40m  15.0  1.78 

20 . 0  1 . 32 

25 . 0  1 . 06 

30 . 0  0 . 90 
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Table    3.6    (continued) 


Saturated  Case 

H 

je 

FS 

5  m 

25 

1  .88 

30 

1  .56 

35 

1  .35 

40 

1  .20 

45 

1  .  12 

50 

1  .04 

55 

0.99 

10  IT. 

15 

1  .95 

20 

1  .47 

25 

1  .  16 

30 

0 .  96 

20  m 

12. 

5 

1  .66 

15. 

0 

1  .37 

17. 

5 

1  .  16 

20. 

0 

0.99 

30  rr. 

10 

1  .79 

12. 

5 

1.41 

15. 

0 

1.15 

17. 

5 

0.97 

40  m 

10. 

0 

1  .61 

12. 

5 

1.28 

15. 

0 

1  .05 
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Interpretation  of  the  Factor  of  Safetu 

This  section  addressee  the  engineering  interpretation 
of  a  computed  factor  of  safety  and  the  choice  of  a  (ninimum 
factor  of  safety  to  achieve  embankment  stability.   The  gen- 
eral rule  of  thumb  in  the  past  has  been  that  a  value  of  the 
factor  of  safety  for  earthuiorks  betueen  1.  3  and  1.  5  is 
acceptable  (Terzaghi  and  Peck<  1967).   The  louei-  limit  is 
for  maximum  loading  conditions  and  the  upper  limit  is  for 
service  conditions  (Meyorhoff.  1970).    If  the  engineer 
designs  an  embankment  slope  with  conventional  methods  of 
analysis  and  these  values  of  the  factor  of  safetyi  it  is 
expected  that  the  embankment  will  perfor-n  satisfactorily. 
Up  to  this  pointi  the  factor  of  safety  that  has  been  dis- 
cussed has  been  a  strength  factor  of  safety.   Since  factors 
such  as  the  dependence  of  the  position  of  the  critical  sur- 
face on  the  factor  of  safety  make  the  interpretation  of  the 
strength  factor  of  safety  difficult*  supplemental  methods  of 
evaluating  the  proximity  of  a  slope  to  lir.it  equilibrium  are 
useful.   Two  such  methods  are    the  geometric  interpretation 
and  the  probabilistic  approach. 

The  Geometric  Interpretation 

The  factor  of  safety  is  usually  defined  as  the  ratio  of 
the  strength  in  the  limit  equilibrium  state  to  the  actual 
stress.   The  choice  of  strength  as  the  variable  for  this 
comparison  is  arbitrary<  however.   It  is  equally  valid  to 
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define  the  factor  of  safety  in  terms  of  some  relevant  dimen- 
sion of  the  slope  geometry.   Consider  a  typical  slope  in 
Figure  3.27a.   Generally^  the  material  parameters  of  the 
soil  J  ♦..  c  and  V  are  known.   The  height.,  H,  and  the  uiidth., 
W>  are  specified  by  need.   The  only  variable  that  is  not 
fixed  is  the  sides  lope.,  B.       It  seems  appropriate,  therefore., 
to  define  the  factor  of  safety  as  the  ratio  of  the  sideslope 
at  limit  equilibrium,  B      ,    to  the  actual  sideslope,  B,    or. 


FS^  =  jB   /je  <3.  11> 

Je     cr 

The  value  of  the  side  slope  at  limit  equilibrium:,  B      ,    may 
be  obtained  graphically  by  plotting  the  strength  factor  of 
saPety  vs.  the  side  slope  for  fixed  values  of  4 ,  c,  >S    and  H. 
j6    is  the  point  on  the  curve  that  corresponds  to  a  strength 
factor  of  safety  of  1.0  (Figure  3.27b). 


Examp le  3. 3 

It  is  desired  to  determine  the  relationship  betuteen  the 
intrinsic  strength  factor  of  safety  of  an  embankment  and  the 

fi      /B   ratio.   The  height  of  the  embankment  is  15.24  m  <50 

cr 

ft).   The  soil  densityi  friction  angle*  and  cohesion  inter- 
cept are  taken  to  be  19.34  kN/m^  CiaS  pcf>,  21°,  and  14.85 
kPd  (310  psf)i  respectively,  to  simulate  long-term  behavior. 
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PROCEDURE  TO  DETERMINE  THE  LIMIT  EQUILIBRIUM 
VALUES  OF  THE  SIDESLOPE  AND  HEIGHT 
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The  first  step  is  to  calculate  the  values  of  the 

intrinsic  strength  factor  of  safety  for  various  side  slopes. 

This  can  be  done  with  the  computer  program  in  Appendix  C. 

The  results  are  given  in  Table  3.7  and  Figure  3.28.   The 

value  of  B         obtained  from  Figure  3.28  is  39.1  .   The  ratio, 
cr 

B      /B,     is  plotted  vs.  the  strength  factor  of  safety  in  Fig- 
ure 3.29.  For  this  case,  the  ratio  B      /B    is  always  greater 

cr 

than  the  value  of  the  strength  factor  of  safety.   One 
shortcoming  of  this  approach  occurs  when  there  is  no  value 
of  sideslope  for  which  the  factor  of  safety  reduces  to  1.0. 
In  such  a  case  the  factor  of  safety  can  not  be  defined  as 
B      /B.       Fortunately,  in  such  cases,  the  shear  strength  is  so 
high  that  a  stability  analysis  usually  is  not  necessary. 

A  factor  of  safety  based  on  the  ratio  of  the  height  at 
limit  equilibrium  to  the  actual  height,  H   /H,  may  be 
defined  in  a  manner  analogous  to  the  development  of  the  side 
slope  factor  of  safety  by  holding  ♦,  c,  V,  and  jB  constant. 
The  actual  geometry  of  the  embankment  and  the  hypothetical 
geometry  at  limit  equilibrium  are    shown  in  Figure  3.30.   The 
graphical  method  for  finding  the  critical  height  is  illus- 
trated in  Figure  3.  28. 


Example  3.  4 

It  is  desired  to  find  the  variation  of  the  ratio  H   /H 

cr 

of  an  embankment  vs.  the  intrinsic  strength  factor  of 
safety.   The  side  slope  of  the  embankment  is  1  to  2 
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Table  3.7   Sideslope  Factor  of  Safety  Calculation 
Examp  1  e  3.  3 


FS 


22 

1.  670 

77  7 

24 

1.  535 

629 

26 

1.  425 

504 

23 

1.  331 

396 

30 

1.  251 

303 

32 

1.  181 

T-,p 
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1.  120 

150 
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1.  069 
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<26.57'^>.   Assume  the  values  of  the  soils  material  parame- 
ters are    the  same  as  in  example  3.3.   The  height  is  vari- 
able.  The  strength  factor  of  safety  uias  obtained  for  vari- 
ous heights  with  the  computer  program  in  Appendix  C.   The 
results  are    given  in  Table  3.8  and  Figure  3.31.   The  criti- 
cal height  is  approximately  44.2  m  (145  ft).   The  critical 
value  is  divided  by  each  height  corresponding  to  a  value  of 
the  strength  factor  of  safety.   The  results  are  given  in 

Table  3.8  and  Figure  3.32.  The  ratio^  H   /H,  is  substan- 

cr 

tially  larger  than  the  value  of  the  strength  factor  of 
safety.   This  is  not  always  the  case.   For  example,  if  IS,    ^, 
and  c  iDere  taken  to  be  19.59  kh4/m^  <124.6  pcf>,  7.3*^,  and 
95.02  kPa  (1984  psf)i  respectively*  to  si<nulate  the  unsa- 
turated short  term  strength  of  St.  Croix  clay  compacted  27. 

uet  of  OMCi  the  results  mould  be  those  shown  in  Figure  3. 33. 

I 
In  this  case  the  difference  between  the  ratio  H   /H  and  the 

cr 

intrinsic  strength  factor  of  safety  is  sxaller.   This  exam- 
ple demonstrates  that  the  strength  factor  of  safety 
correspoTid'^  to  different  values  of  the  H   /H  ratio,  depend- 
ing on  the  values  of  ♦  and  c.  An  analagojs  remark  applies  to 
the  side  slope  ratio. 


The  Probabilistic  Approach 

The  conventional  definition  of  the  factor  of  safety 
overlooks  the  variability  of  material  parameters  because  it 
uses  deterministic  (single-valued)  soil  properties.   The 
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Table  3.8   Factor  of  Safety  Calculated  Based  on  Height 
Criterion  -  Example  3.4 


H  <ft> 

H,,/H 

FS 

10 

14. 

500 

3.  456 

20 

7. 

250 

2.  192 

25 

5. 

800 

1.  932 

30 

4. 

833 

1.  759 

35 

4. 

143 

1.  63 

40 

3. 

625 

1.  533 

45 

3. 

222 

1.  458 

50 

2. 

90 

1.  393 

55 

2. 

636 

1.  347 

60 

2. 

417 

1.  304 

65 

2. 

231 

1.  268 

70 

2. 

071 

1.  235 

75 

933 

1.  207 

80 

813 

1.  183 

85 

706 

1.  166 

90 

611 

1.  143 

95 

526 

1.  125 

100 

450 

1.  113 

105 

380 

1.  099 

110 

318 

1.  086 

115 

261 

1.  070 

120 

208 

1.  060 

125 

160 

1.  050 

130 

115 

1.  035 

135 

074 

1.  025 

140 

036 

1.  010 

105 


FS 

i, 
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FIGURE  3.3 1       DETERMINATION  OF  Hjj^    -    EXAMPLE  3.4 
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FIGURE   3.32      FACTOR  OF  SAFETY  BASED  ON  THE  HEIGHT 
CRITERION  VS.  THE  STRENGTH  FACTOR  OF 
SAFETY -EXAMPLE  3.4 
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variability  is  only  taken  into  account  through  the  personal 
judgment  exercised  in  the  selection  of  the  soil  properties. 
Variability  in  soil  properties  arises  through: 

(1)  material  and  sample  non-homogeneities  that 
do  not  represent  the  whole  soil  ("true" 
variability ) 

(2)  sampling  errors  caused  by  disturbance 
during  the  sampling  process 

(3)  errors  that  occur  when  tests  are  not 
performed  according  to  a  standard. 

Uncertainty,  however,  is  not  limited  to  the  variability 
observed  in  the  basic  variables.   Analytical  models  and 
laboratory  and  field  experiments  are    often  only  an  idealized 
representation  of  reality.   Predictions  made  on  the  basis  of 
these  models  and  experiments  may  be  inaccurate  and  contain 
uncertainty.   Therefore,  the  capacity  of  a  slope  to  resist 
loading  will  not  have  a  unique  value.   Similarly,  the  load 
(or  demand)  on  the  trial  surface  will  have  a  distribution  of 
values.   These  distributions  can  be  represented  by  probabil- 
ity density  functions.   The  probability  that  a  slope  will 
reach  a  state  of  limit  equilibrium  equals  the  probability 
that  the  capacity  will  be  less  than  tite  demand  (Yao,  1982), 
i.  e.  , 

S 
max 

P-.   =   p<:R<S>    =      X        f--<s>      F_<s>    ds  <3.12> 

^  R    .         ^  "^ 

mm 
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where 

R  denotes  strength  (or  capacity) 

S  denotes  load  < or  demand) 

f_  is  the  probability  density  function  <PDr>  of  the  load 

F_  =  cumulative  distribution  function  <CDr>  of  the 

resistance. 

R  .   =  minimum  value  of  the  strength  <see  Figure  3.34> 
mm 

S     =  maximum  value  of  the  load  <see  Figure  3.34> 
max 


Capacity-demand  problems  are    simplified  in  the  case  of 
slope  stability  analysis  because  the  capacity-demand  ratio 
for  a  slope  of  specified  geometry  and  material  parameters  is 
identical  to  the  factor  of  safety.   Therefore,  the  probabil- 
ity of  failure^  P.,  is  also  equal  to  the  probability  that 
the  factor  of  safety  is  less  than  one.  oc 

P^  =  PCFS  <  1 .0>  <3. 13> 

This  calculation  is  performed  by  integrating  the  probability 
density  function  of  the  factor  of  safety  up  to  a  value  of 
one.  The  density  function  of  the  factor  of  safety  depends  on 
tha  distributions  of  the  soil  density  and  shear  strength 
variables.   When  the  distribution  of  these  variables  is 
knou/n.  the  density  function  of  the  factor  of  safety  may  be 
obtained  by  Simulation  (Yao.  1982). 

Simulation  is  essentially  a  controlled  statistical  sam- 
pling technique  which  can  be  used  to  study  complex 
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S  min 


R  min      "R  S  max     R  max 


FIGURE    3.34       PROBABILITY  DENSITY  (OR  DISTRIBUTION) 
FUNCTIONS  OF  CAPACITY  AND  DEMAND 


Ill 

stochastic  systems  when  analytical  and/or  numerical  tech- 
niques do  not  suffice.   A  necessary  part  of  any  such  pro- 
cedure is  an  algorithm  for  random  number  generation.   A  ran- 
dom number  generator  produces  sequences  with  probability 
density  functions  that  are    uniformly  distributed  between  0 
and  1  and  that  possess  the  appearance  of  randomness.   Most 
computer  systems  have  a  built-in  random  number  generator. 
The  inverse  transformation  method*  often  called  the  Monte- 
Carlo  method  of  simulation*  is  used  to  generate  non-uniform 
random  number  X»  with  cumulative  distribution  function 
F  <x> .   The  algorithm  is  very  simple: 

1.  Generate  number  Q  uniformly  distributed  between 
O  and  1. 

2.  Return  X  =  f"^<:Q> 

X 

<F         is  the  inverse  function 

X 

corresponding  to  F  > 
This  alqorithm  assumes  that  the  equation  F  <x>  =  Q  can  be 

X 

solved  explicitly.   Other  distributions  can  be  simulated  by 
direct  generation  methods  such  as  composition  methods  and 
rejection-acceptance  methods  (Fishman.  1978).   The  probabil- 
ity of  failure  equals  the  cumulative  distribution  function 
of  the  factor  of  safety  up  to  a  value  of  one>  i. e. * 


P^  =  F^_  CFS  =  1.0>  <3.14> 

r     rs 

Since  the  Monte-Carlo  method  uses  the  actual  distributions 
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of  the  variables  that  affect  the  factor  of  safetg.  it  can. 
in  principle,  generate  the  actual  density  function  of  the 
factor  of  safety.   Unfortunately,  this  requires  a  very  large 
number  of  simulations  and  consequently,  a  great  amount  of 
computational  effort. 

A  simpler  approach  used  in  practice  to  obtain  the  pro- 
bability of  failure  is  to  assume  a  distribution  of  the  fac- 
tor of  safety  and  to    generate  statistical  moments  of  this 
distribution  from  the  statistical  moments  of  the  dependent 
variables  (density  and  shear  strength  parameters).   Two  pro- 
cedures to  obtain  these  moments  are  the  Taylor  Series  expan- 
sion method  and  the  Point  Estimates  Method  <Yao,  1982). 

The  Taylor  Series  expansion  requires  partial  deriva- 
tives of  the  factor  of  safety  with  respect  to  each  of  the 
variables  that  affect  the  factor  of  safety  (Harr.  1977). 
The  difficulty  of  evaluating  these  derivatives  limits  the 
usefulness  of  the  Taylor  Series  expansion  for  modelling  the 
density  function  of  the  factor  of  safety. 

The  Point-Estimates  Method  is  an  approximate  procedure 
for  calculating  the  statistical  moments  of  the  factor  of 
safety  which  does  not  require  evaluation  of  the  derivatives 
of  the  factor  of  safety  relative  to  the  variables  affecting 
it  (Rosenblueth.  1975).   This  method  approximates  the  sta- 
tistical moments  of  a  function  Y  (i.e..  the  factor  of 
safety)  by  "replacing"  the  distribution  of  the  variables 


i.i;i 


a-f'fvctlng  th»  function  with  point  •ttint«tcfi  (or  uitights)  at 

properly  selected  values  of  the  variables.   For  example,  the 

th 
estimated  value  of  the  n    statistical  moment  of  the  func- 
tion Y  that  depends  on  two  random  variables  may  be  feund 
with  the  following  two-point  procedure: 


EO'''^  =  p++y++  +  p+_yr-  *  P-+y-+  *  P— y-  <3.i5> 


y^i  =  YCX^  ±c^j,  X^    ±a^^l  <3.16> 

and 


P  =  weighting  factor 
X..,X^  =  random  variables 
*'xi'''^X2  ~  standard  deviations  of  the  random  variables 
X  =  mean  value  of  the  variable/  X 


If  X   and  X^  have  symmetric  distribut  ion=:^  the  lieighting 
factor  is  1/4.   The  point  estimates  method  can  be  adapted  to 
handle  any  number  of  correlated  and/or  uncorrelated  random 
variables  (Rosenb lueth/  1975). 

To  illustrate  the  application  of  the  probabilistic 
approach  to  slope  stability  problems<  the  friction-circle 
slope  stability  program  (Appendix  C)  was  adapted  to  calcu- 
late the  probability  of  failure  of  a  sinple  slope.   The 
resulting  program  is  given  in  Appendix  E.   This  program 
assumes  that  the  material  parameters  defining  the  soil 
stf'ength  and  density  have  symmetric  distributions.   These 
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dietributions  are   defined  with  mean  values,    coefficients  of 
variation^  and  upper  and  lower  bounds.   The  material  parame- 
ters are   assumed  to  be  uncorrelated.   The  probability  den- 
sity function  of  the  factor  of  safety  is  assumed  to  be  beta 
distributed  (Appendix  D).   The  mean  and  variance  <i.e. >  the 
first  two  statistical  moments)  of  the  factor  of  safety  are 
obtained  with  a  two— point  point  estimates  procedure.   The 
lower  bound  of  the  distribution  of  the  factor  of  safety  is 
the  factor  of  safety  that  is  obtained  using  the  maximum 
value  of  the  soil  density  and  the  minimum  values  of  the 
Mohr-Coulomb  parameters ,  ♦  and  c .   The  upper  bound  is  the 
factor  of  safety  that  is  obtained  using  the  minimum  value  of 
the  density  and  maximum  values  of  the  Mohi — Coulomb  parame- 
ters.  Finally*  the  probability  of  failure  is  calculated  by 
integrating  equation  3.  13  numerically. 

Example  3.  5 

It  is  desired  to  determine  the  probability  of  failure 
of  the  slope  shown  in  Figure  3.35.  The  average  ratio  VH/c 
is  25/3 

where 

V  =  the  mean  value  of  the  soil  density 

c  ■=  the  mean  value  of  the  cohesion  intercept 

of  the  soil 
H  =  the  height  of  the  slope 

The  variability  of  the  material  parameters  is  given  in  Table 
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3.  9.   The  mean  factor  of  safety  corresponding  to  the  mean 
values  of  the  soil  density  and  the  Mohi — Coulomb  strength 
parameters  is  found  to  be  1.28  using  the  friction  circle 
method.   The  program  in  Appendix  E  is  used  to  quantify  the 
effect  of  the  variability  of  the  cohesion  intercept  on  the 
probability  of  failure.   The  results  are   reported  in  Table 
3.  10  and  Figure  3.  36.   An  increase-  in  the  variability  of  the 
cohesion  intercept  increases  the  probability  of  slope 
failure.  Therefore*  slopes  uith  equal  mean  factors  of  safety 
are  not  necessarily  equally  safe.  Historical  records  indi- 
cate that  earth  dams  designed  with  ordinary  techniques  for  a 
factor  of  safety  of  1.3  to  1.5  have  a  rate  of  failure  on  the 
order  of  one  in  two  thousand  (Meyerhof/  1970).   However> 
further  work  needs  to  be  done  in  order  to  recommend  design 
values  of  the  probability  of  failure  for  various  situations 
such  as  sudden  drawdown  and  long-term  conditions. 

Frequently*  the  probability  of  failure  that  is  calcu- 
lated based  on  ordinary  estimates  of  material  variability  is 
much  higher  than  the  values  reported  by  Meyerhoff.   This 
discrepancy  arises  because  designers  typically  use  a  lower 
bound  on  their  strength  parameters  rather  than  mean  values 
when  they  calculate  the  factor  of  safety.  If  mean  values  of 
the  strength  parameters  were  used*  a  higher  factor  of  safety 
has  to  be  employed  to  assure  adequate  performance. 
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Table  3.9   Material  Parameters  -  Example  3. S 


Parameter 


c 

min 
c 

max 


Value 


5'=' 

0 

.6   ♦ 

1 

.4   4 

B5    V 

1 

.15   y 

e>6   c 

1 

.33   c 

Parameter    Coef f ic lent  of  Variation 


0.2 

0.1 

variab le 
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Table  3.  10   Probability  of  Failure  vs.  Variability  of  Cohesion  - 
Examp  le  3.  5 


FS 


20°  1 . 23 


1  .05 


^^c 

Pf 

.01 

.016 

.  05 

.  023 

.  10 

.  047 

.  20 

.  145 

.  30 

.  253 

.  33 

.  284 

.01 

.  226 

.  05 

.  292 

.  10 

.  330 

.  20 

.  401 

.  30 

.  449 

.  33 

.  462 

f>   =   si  des  1  ope 

FS  =  f-actor  of"  safety  based  on  the  mean  value  of  +..  Cj    and 

u   =  coefficient  of  variation  of  the  cohesion 

p"  =  probability  of  failure. 
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FIGURE  3.36       PROBABILITY  OF  FAILURE  VS.  COEFFICIENT 
OF  VARIATION  OF  THE  COHESION 
INTERCEPT  -  EXAMPLE  3.5 
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IV  -  SETTLEMENT  CONSIDERATIONS  FOR  EMBANKMENT  DESIGN 

Once  an  embankment  has  been  checked  for  intrinsic  and 
overall  stability  of  its  sideslopes>  the  settlement  of  the 
embankment  should  be  investigated.   Embankment  settlement  is 
comprised  of  displacements  that  occur  within  the  embankment 
itself  as  well  as  compression  of  fine  grained  soil  lagers 
that  underlie  the  embankment.   Displacements  within  the 
embankment  itself  are  caused  principally  by: 

1)  partially  saturated  compression  under  the  body  forces 
of  the  fill.   This  occurs  as  rapidly  as  the  fill  is 
constructed  (DiBernardo>  1979). 

2)  volume  change  due  to  an  increase  in  moisture  content 
of  the  compacted  embankment.   This  is  thought  to  be 
the  major  source  of  displacements  in  compacted  clay 
embankments. 

Compression  of  the  fine  grained  soil  layers  below  the 
embankment  consists  of: 

1)  immediate  settlement  that  occurs  at  constant  volume. 
Immediate  settlement  is  completed  at  the  end  of 
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embankment  conetruction  (consequently*  it  has  no 
effect  on  the  embankment  performance  and  may  be 
neglected). 

2)  time  dependent  consolidation  settlement  that  occurs 
as  excess  pore  pressure  induced  by  the  embankment 
construction  dissipates 

3)  secondary  compression  settlement  that-  occurs  after 
time  dependent  consolidation  is  complete  (calcula- 
tion of  the  secondary  compression  settlement  is 
frequently  omitted  because  it  is  assumed  to  be 
small  compared  to  the  consolidation  settlement). 

Saturation  Induced  Displacement  of  Compacted  Clau  Embankments 

When  a  sample  of  compacted  clay  becomes  saturated*  it 
may  either  swell  or  settle  depending  on  the  mineralogy  of 
the  clay>  the  presaturation  water  content*  the  compactive 
effort*  and  the  overburden  pressure.   Prediction  models  have 
been  developed  for  the  volume  change  due  to  saturation  for 
laboratory  compacted  St.  Croix  clay  (DiBernardo*  1979).   An 
extension  of  this  model  for  field  compacted  soils  with  a 
range  of  plasticity  index  values  was  proposed  by  Lin  (1981). 
However*  since  the  data  base  for  these  models  is  limited*  it 
is  recommended  that  the  volume  change  due  to  saturation  be 
estimated  from  tests  on  the  soil  that  will  be  used  in  the 
embankment. 
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To  insure  that  the  testing  will  best  simulate  the 
volume  change  in  the  embankment  due  to  saturation,  the  fol- 
lowing guidelines  should  be  followed: 

1)  The  soil  sample  should  be  compacted  to  the 
expected  state  of  compaction  in  the  field  according 
to  the  procedure  described  in  Chapter  II. 

2)  The  presatiiration  water  content  of  the  soil  should 
be  the  same  as  is  expected  in  the  actual  embankment. 

3)  The  soil  sample  should  be  subjected  to  a  pressure 
equal  to  the  overburden  pressure  in  the  embankment. 
Therefore,  to  simulate  the  variation  of  volume  change 
with  depth,  the  test  samples  must  be  subjected  to 

a  range  of  overburden  pressures. 

4)  Measurement  of  volume  change  of  the  soil  sample  should 
be  made  from  the  time  that  the  soil  is  back-pressure 
saturated  until  the  volume  change  ceases.   For  de- 
tails of  the  test  procedure,  see  DiBernardo  (1979). 

If  it  is  assumed  that  all  volume  change  occurs  verti- 
cally, it  is  possible  to  estimate  the  settlement  of  the  sur- 
face of  the  embankment  with  the  following  expression: 


S  =  T^  S   U<z>  dz  <4.1a> 


where 


S  -  the  settlement  due  to  saturation 
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U<z)  =  the  '/.   volume  change  at  the  depth  2 

z  =  the  depth  beneath  the  top  of  the  embankment 
H  =  the  height  of  the  embankment 

If  the  distribution  U(z)  can  be  idealized  as  a  sequence  of 
strata  each  uiith  its  own  uniform  value  of  U,  then  equation 
4. la  may  be  evaluated  numerically  with  the  following  expres- 
sion 


1    ^ 
S  =  7^   E   U.  Az.  <4.1b> 

1  =  1 


where 


U.  =  the  v.   volume  change  of  stratum  i 

6z .  =  the  thickness  of  stratum  i 

1 

n  =  the  number  of  strata 


It  is  interesting  to  note  that  it  is  possible  for  the  upper 
portion  of  an  embankment  to  be  swelling  while  the  lower  por- 
tion is  settling.   In  fact*  it  is  theoretically  possible  to 
build  an  embankment  whose  net  saturation  settlement  is  nill 
by  specifying  the  compaction  so  that  there  are  compensating 
zones  of  swelling  and  settling  soil  within  the  embankment. 
6enerally>  however/  it  is  best  to  design  the  embankment  to 
settle  slightly  because  swelling  of  soil  beneath  the  road 
bed  can  cause  severe  pavement  distress. 
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Consolidation  Settlement  of  Compressible  Soil 
Lauers  Beneath  the  Embankt^gnt 


Magnitude  of  Settlement 

When  a  saturated  clay  sample  is  subjected  to  an  axial 
stress  change  in  a  standard  consolidation  test,  an  excess 
pore  pT-essure  equal  to  the  stress  change  is  induced  in  the 
sample.   As  time  proceeds,  the  excess  pore  pressure  uill 
dissipate  and  the  sample  luill  settle  by  a  volume  equal  to 
the  volume  of  the  dissipated  pore  water.   The  relationship 
between  the  axial  stress  and  the  void  ratio  after  all  of  the 
excess  pore  pressure  has  been  dissipated  is  described  in 
Figure  4.  la.   For  purposes  of  analysis,  the  relationship  may 
be  replaced  uith  the  representation  shown  in  F'igure  4.  lb 
U'hich  1=    de-fined  by  the  compress  iori  index,  C.  ..  the 
recoiTipTe-ssion  index.,  C  .,  and  the  preconsol  id-at  ion  pressure. 


If  the  exist  iriQ  axial  pressure,  cf'  ,  is  equal  to  the 
pr ec onso 1 1  da t i on  pressure,  the  soil  is  normally  consoli- 
dated.  The  settlement  of  the  soil  due  to  the  axial  pressure 
change  is: 


1  +  e 


C  '  loci 


a  •'   +  Act 

V 


•;;  4  . 2 : 


where 

S  =  the  settlement  of  the  clay  layer 
H  =  the  thickness  of  the  clay  layer 


125 


CO 


o 

> 


»   log  o 
Logarithm  of  Axial  Stress 
FIGURE  4.1a  TYPICAL  e-log  a'  RELATIONSHIP 


cd 
IT 


o 
> 


-^  log  C 


Logarithm  of  Axial  Stress 

FIGURE  4. 1  b  SIMPLIFICATION  OF  TYPICAL  e-log  o' 

RELATIONSHIP  FOR  ANALYTICAL 
PURPOSES 
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dj<7    -    the  increase  iri  axial  stress  of  the-  'lay  len^ier 
e   =  the  initial  void  ratio  of  the  c 1 au  1 auer 
a'   =  +  he  existing  ov^erburden  pressure 

If  the  existing  overburden  pressure  i&  less  than  the 
preconsolidation  pressure,  the  soil  is  overconsolidated 
Ihe  settlement  due  to  the  axial  pressure  change  is; 


-     -  CC   loq  ;;: 2 

1  +  e     r    -    c 


<:  4  . ::;  ::■ 


liihen    <c         +    iicr>    <    o"'     ,     or 
V  -        p 


<j  '       "•"    li  I?    —    c ' 


— CC    •  loo    — 

1    +    e  r         ~    <j 


C:       loci 


■3        '::4.4: 


inhon    Cc         +    in-;' 

V 


If  the  existing  overburden  pressure  is  greater  than 
'^"'/    ^^'''='    soil  is  underconsol  i  dated  .   The  set  t  1  e;T.erit  due  to 
the  axial  pressure  change  is: 


1  +  e 


C   lo.: 


':;4 


These  expressions  for  settlement  are  exact  provided 


that 


1  '  the  values  of  e  .  C  ,  C 
o'   r    f. 


O"  ,.  and  a    are  the 
V         p 


same  in  the  consolidating  layer  as  measured  in 
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the  consolidation  test 

2)  the  stress  change  in  the  consolidating  layer  does 
not  vary  with  depth 

3)  the  excess  pore  pressure  caused  by  the  embankment 
loading  is  the  same  as  in  the  consolidation  test. 

Deviations  from  the  first  two  assumptions  may  be  accounted 
for  by  dividing  the  consolidating  layer  into  artificial 
strata.   Deviations  from  the  third  assumption  may  be 
accounted  for  by  use  of  a  correction  factor  that  will  be 
discussed  later. 

To  determine  the  value  of  e   in  each  stratum,  one  can 

o 

assume  that  the  difference  in  the  initial  void  ratio  between 
the  center  of  the  stratum  and  the  position  of  the  soil  sam- 
ple, Ae  ,  will  be  equal  to  the  change  in  void  ratio  that 
would  occur  if  the  soil  sample  were  to  move  along  the  e  - 
log  ct'  curve  by  a  stress  change  equal  to  the  difference  in 
overburden  pressure  between  the  position  of  the  soil  sample 
and  the  center  of  the  stratum. 

If  the  center  of  the  stratum  is  above  the  elevation  of 
the  sample,  6e   will  be  positive.   If,  in  addition,  the  soil 
is  overconsol i dated  above  the  elevation  of  the  sample 


Ae   =  C   log  Ca     ,  /a'   >  <:4.6a> 

o     r        sample    vo 


where 
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a     ,   =  the  overburden  pressure  in  the  soil  sample 
sample 

a'    =  the  overburden  pressure  in  the  stratum 
vo 


If  the  soil  above  the  elevation  of  the  san:ple  is  considered 
to  be  underconso 1 i dated : 


ie   =  0  <4.6b> 

o 

If  the  soil  is  underconsol idated  at  pressures  above  the 

preconsol  idat  ion  pressure  of  the  sample >  o"'  ,    and  overconso- 

1  idated  at  pressures  be  Ion*  C 

P 


iie=Clo9<cr'/a'       >  <4.6c> 

o  r         ^  p        vo 


If    the    soil    is    normal lu   consolidated   at    pressures    above    o' 

P 

and  overconsol idated  at  pressures  beneath  o" 


£ie   =  C   log  Ctj'  ,  /cr'  >  +  C   log  C^"    /o'       >        <4.6d> 

o     c    ^     sample    p      r    ^     p    vo 


If  the  center  of  the  stratum  is  below  the  elevation  of 

the  sample^  ie   u^ill  be  negative.   If  the  soil  is  overconso- 
o 

lidated  beneath  the  elevation  of  the  sasr.ple 


Ae   =  -  C   log  Ccr'   /cr'     ,  >  <4.fce> 

o      r    ^     vo    sa-i^ple 


If  the  soil  is  overconsol  idated  at  pressui-es  below  q"      and 

P 

normally  consolidated  at  pressures  above  a'    , 

P 
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Ac      =    —    r     '  IriQ    im'    /a  ,     >    —    C     '  ]  oq    ^a'       /'T'     > 

"o  "r  -  p       satTiple  c  -  vo         pi 


<4  .  t.f 


If  the  soil  is  overc onso 1 i dated  throughout  its  stress  path 


C  ■  loQ  <.<7'    /(J  .     ':• 

r    ~     p      saiTiplt- 


'-.4  .  tq 


If  the  soil  is  underc onsol idated  throughout  its  stress  path 


;  4  .  6h-i 


Like    the    void    ratio,     a'       '•••aries    with    depth    in    a    c  1  au 

P 

layer-.   Typa  cal  ly  ^  ct'   decreases  aiith  deptPi  i-.mt  i  1  the  sc'il 

is  normally  consolidated.   Thereafter  it  assu^-nes  a  value 

e:iual  to  the  overburden  pressure  (Holtz  and  Kovacs/   19S1). 

This  15-  illustrated  in  Fiqure  4.2a.   The  variation  of  a 

P 

with     depth    may     be    represented     by     tujo    straight -line     segaients 

''FiQi.ire    4    I'P' .        ijnce    the    values    of    <-r  at    tl-ie    endp'Oirit^     r.f 

P 

th-ie    tii'o    seoiTients    are    specified,     the    valLie    of    o''       iTiag    be 

P 

obtained  by  interpolation  in  any  stratum.    Developing  the 

o'       ciTof  lie  reqi-iires  more  consolidation  te~''s  thati  are  nor- 
P 

(Tially  run.   This  is  not  prohibitive-  hD'jjever-  because  these 
consolidation  tests  can  be  I'un  as  constant  gradient  tests 
and  completed  in  one  tuork  day. 


The  change  in  vertical  pressure  due  to  the  embankment 

load  can  be  approximated  with  the  expressions  presented  in 

Hppendi.x  F.   C_  and  C   are  assumed  constant  for  a  giv^'eri  clay 
layer. 
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Overburden  Pressure 
Profile 


Preconsolidation  Pressure    Profile 


FIGURE  4.2a 

TYPICAL  PRECONSOLIDATION  PRESSURE  PROFILE 


Sand 


Simplified  Preconstruction 
Pressure   Profile 

Clay 


FIGURE  4.2b 

SIMPLIFIED  PRECONSOLIDATION  PRESSURE  PROFILE 
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Once    the   values    ofe,C,C,a'',a''       and   Aa    are 

ore     V     p 

determined  in  each  stratum,  the  settlement  of  the  consoli- 
dating layer  may  be  calculated  by  evaluating  equations  4. 2  - 
4.5  for  each  stratum.   The  computer  program  in  Appendix  G 
has  been  provided  to  facilitate  computations. 

Examp  le  4.  1 

It  is  desired  to  make  a  preliminary  estimate  of  the 
consolidation  settlement  of  the  embankment  built  over  the 
clay  layer  shouin  in  Figure  4.3  without  performing  consolida- 
tion tests.  The  compression  index  may  be  estimated  using  an 
appropriate  correlation  with  standard  index  tests  (Terzaghi 
and  Peck.  1967).   The  recompression  index  is  assumed  to  be 
one  tenth  of  the  compression  index.  For  purposes  of  this 
example,  ct'   was  taken  to  be  95.76  kPa  <2000  psf)  at  the 
depth  of  the  sample.   cr'   was  assumed  to  equal  the  overbur- 
den pressure  when  the  overburden  pressure  exceeds  95.76  kPa. 

It  is  instructive  to  consider  the  stress  changes  in  the 
clay  layer  caused  by  the  embankment  before  actually  calcu- 
lating settlement.   This  was  done  with  the  computer  program 
in  Appendix  F.   The  results,  which  are  shown  in  Table  4. 1. 
indicate  that: 

1)  The  increase  in  vertical  stress  at  a  given  depth 
in  the  clay  layer  is  approximately  constant 
across  the  embankment 

2)  The  vertical  stress  increase  does  not  attenuate 
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|- 2  1. 34m  (70ft) — *| 


2. 1  9m  (40ft)    ^1     j)<=,9.64kN/m3(,25pcf^ 
2  Compacted  Clay 


1.52m  (5ft) 


Sand 


1.52m  (5ft) 


5.24m  (50ft)      Clay 


"77777 77777- 


4  =l7.28kN/m3(l  lOpcf) 

^     

Cr  '  8-07kN/m3  (115  pcf) 


^  =l9.64kN/m3(i25  pcf) 

sat  ^ 

LL=45 
Cc   =0.32 
Cr   =0.03 

Op  =  95.76kPa  (2000  psf) 
fio  =0.7 


"77777^— 777py- 


FIGURE  4.3 


EMBANKMENT  -  EXAMPLE  4. 


Table  4.1   Stress  Change  Beneath  Embankment 
Figure  4.  3 
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X  <:  m  > 


0. 

00 

3. 

05 

0. 

00 

4. 

57 

0. 

00 

6. 

10 

0. 

00 

7. 

62 

0. 

00 

9. 

14 

0. 

00 

10. 

67 

0. 

00 

12. 

19 

0. 

00 

13. 

72 

0. 

00 

15. 

24 

0. 

00 

16. 

76 

0. 

00 

18. 

29 

10. 

67 

3. 

05 

10. 

67 

4. 

57 

10. 

67 

6. 

10 

10. 

67 

7 

62 

10. 

67 

9. 

14 

10. 

67 

10. 

67 

10. 

67 

12. 

19 

10. 

67 

13. 

72 

10. 

67 

15. 

24 

10. 

67 

16. 

76 

10. 

67 

18. 

29 

(u(j    CkPa. 


239.  0 

238.  0 

236.  3 

233.  8 

230.  7 

228 

222 

218 

213 

208 

203 


229.  9 
225.  1 
220.  3 
215.  6 
210.  9 
206.  1 
201.  4 
196.  8 
192.  2 
187.  7 
183.  2 


X 

At   <kPa:.' 
xz 

173.  7 

-74.  4 

163.  4 

-72.  3 

153.  8 

-69.  6 

69.  0 

-66.  6 

61.  2 

-63.  2 

54.  3 

-59.  9 

48.  3 

-56.  5 

42.  9 

-53.  2 

38.  2 

-50.  0 

34.  0 

-47.  0 

30.  4 

-44.  1 

148.  8 

8.  0 

118.  5 

8.  0 

94.  4 

6.  5 

75.  4 

4.  5 

60.  6 

2.  3 

49.  1 

0.  5 

40.  1 

-1.  0 

33.  0 

-2.  1 

27.  4 

-2.  9 

22.  9 

-3.  5 

19.  3 

-3.  9 
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appreciably  with  depth. 

3)  The  increase  in  horizontal  stress  is  small  com- 
pared to  the  increase  in  vertical  stress  every- 
where except  near  the  top  of  the  clay  layer. 

4)  Although  the  clay  layer  is  moderately  overcon- 
solidated.  the  vertical  stress  increase  occurs 
principly  at  values  higher  than  the  preconsolida- 
tion  pressure. 

Since  most  of  the  parameters  affecting  the  consolida- 
tion vary  with  depth  in  the  clay  layer,  the  layer  must  be 
divided   into  a  sufficient  number  of  strata  to  insure  ade- 
quate accuracy.   This  can  be  achieved  with  the  computer  pro- 
gram in  Appendix  G.   The  settlement  at  the  centerline  of  the 
embankment  shown  in  Figure  4.3  is  given  in  Table  4. 2  as  a 
function  of  the  number  of  strata,  n.   Five  or  more  strata 
will  suffice  in  this  case.   The  settlement  at  the  embankment 
centerline  is  1.  3  m.   Once  the  number  of  strata  to  be  used 
in  the  analysis  is  known,  the  lateral  variation  of  the  set- 
tlement across  the  embankment  can  be  calculated.   The 
results  are    shown  in  Table  4.3. 

As  noted  before,  the  difference  in  the  stress  paths  of 
the  consolidation  test  and  the  clay  layer  beneath  the 
embankment  will  give  rise  to  the  generation  of  different 
amounts  of  excess  pore  pressure,  and  consequently  to  dif- 
ferent amounts  of  settlement.   An  approximate  method  for 
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Table  4.2   Settlement  vs.  Number  of  Strata  -  Example  4.  1 


S(m) 


1.  372 
1.  305 
1.  295 
1.  308 
1.  295 
1.  298 
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Table  4.3   Settlement  and  Differential  Settlement  Along 
Profile  of  Embankment  -  Example  4.  1 


x(m)         S(m) 


0. 

000 

298 

1. 

524 

298 

3. 

048 

292 

4. 

572 

289 

6. 

096 

280 

7. 

620 

268 

9. 

144 

250 

10. 

668 

228 
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dealing  with  this  discrepancy  is  discussed  m  the  follouiru 
paragrap  h  s. 

The  settlement  of  a  compressible  layer  can  be  defined 


■field 


a  ■  Li '  dz 


:  4  .  S 


where 

a^   =  th"?  coef'f  ic  ierit  of  vertical  coaipressibi  1  i  ty  of  the 
soil 

u  =  the  excess  pore  pressure  due  to  loading 

H  =  the  thickness  of  the  soil  layer 
dz  =  the  thickness  increment  in  the  soil  layer 

In  the  f  1 1- 1  d  the  excess  pressure  will  be  (  H  o  ]  t  7  and  ^  o  v  a  c  s  , 
19B1): 


u  =  EC licr _  _   +  a •  ir 


:  4  . 9 ; 


wh  ere 


':.■  O        +    rr„  +  rr   )  /3 
1      ?      3 


=  T  C  '■  '7 


+  '■.  "7  _  —  a  _  >   +  '■  iT . 


'^l'     ' 


1/i 


=  principal  stresses: 
a-E  =  pore  pressure  parameters  that  are    determined 
experimentally. 
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The  excess  pore  pressure  generated  in  a  consolidation 
test  is  equal  to  the  change  in  vertical  stress,  i. e. . 


u  =  <ia  <4.10> 

V 


It  follouis  that  the  value  of  the  settlement  in  the 
field  will  be 

H 

S^.  ,^  =  /  a  *B  Cia   .  +  air   ^3  dz       <4.11> 
field    ^   V      oct        oct 

and  that  the  value  of  settlement  if  excess  pore  pressures 
are    equal  to  those  in  a  consolidation  test  will  be 

H 
S,  ^  =  X  a  '^^      dr  <4.12> 

lab    o   V    ^ 

By  taking  the  ratio  of  the  settlements  in  equations  4.11  and 
4.12/  a  correction  factor  to  apply  to  settlements  computed 
with  the  laboratory  pore  pressures  may  be  obtained  (Skempton 
and  Bjerrum*  1957).   This  correction  factor  is: 


H 

X   -a      BCia      .    +   a"6T      .3    dz 

Afield        O      ^  ^^^  °^*  ,.    .3, 

^   =    — <4.13> 

Jj       V  V 


Ordinarily<  the  parameteri  B>  is  1.0  for  saturated 
clays.   The  parameter^  a>  depends  on  the  stress  path.  The 
results  of  Example  4.  1  show  that  the  stress  path  beneath  the 


I3<» 


tmb«nl(m»nt  it  •«t*ntlallu  trlaxial  comprtculon    For  thl« 
case>  the  parameteri  a,     may  be  obtained  uith  the  following 
expression 

a  =  3<A    -  i>/\l2  <4.14> 

ac    3 

where  A    =  iu/<i"7<  -  <!la^>,  is  measured  in  a  tri axial  test, 
ac  1      o 

The  parameteri  a<  depends  on  the  amount  of  strain.   When  the 
soil  is  treated  as  an  isotropic  elastic  solid  (which 
incidentally/  is  the  assumption  that  is  made  in  calculating 
the  stress  changes  caused  by  the  embank-nent )  the  parameter/ 
ai     is  zero. 

Assuming  that  the  stress-strain  behavior  of  a  soil  may 
be  idealized  as  in  Figure  4. lb/  the  value  of  the  coefficient 
of  axial  compressibility  in  the  normally  consolidated  range 
fDay  be  taken  to  be 


C  /In  10 
^v  =  CI  '  e>  a-  '■''■'-•' 


Equation  4.13  may  be  simplified  by  assuming  that  a  is 
constant  and  that  B  equals  unity.  The  resulting  expression 
is 


H 

S    <i<T   .  +  air   .  >  dz 
i.  oc  t        oc  t 

M  = ^ <4.16> 

X  icr   dz 
0         ^ 
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Assuming  that  the  soil  may  be  considered  to  be  an  isotropic 
elastic  solid/  this  may  be  further  simplified  to 


H 


t*   = 


2   /    <A<7j    +   iOg   +   ia2>dz 


0 


(.4.  17> 


Examp le    4. 2 

It  is  desired  to  calculate  the  settlement  correction 
factor  for  the  centerline  of  the  embankment  in  example  4. 1. 
The  values  of  the  stress  changes  due  to  the  embankment  load 
(obtained  u/ith  the  computer  program  in  Appendix  F)  are 
presented  in  Table  4.4.   Using  equation  4.  17»  the  resulting 
correction  factor  is  0.66.   This  means  that  the  expected 
settlement  uiill  be: 


Afield  =  ^^'^lab  =  ^  ^^   ^  ^■^^■'•'   '      Q^" 


At  positions  other  than  the  centerline.  the  stress  changes 
and  hence/  the  correction  factor/  will  be  slightly  dif- 
ferent. 

Time-Rate  of  Settlement 


When  the  magnitude  of  consolidation  settlement  is  large 
enough  to  be  of  concern,  it  is  uiorthwhile  to  predict  how 
much  of  this  settlement  will  occur  during  the  service  life 
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of  the  embankment.   For  a  homogeneous  soil  layer  beneath  a 
long  linear  embankment,  consolidation  occurs  due  to  the  dis- 
sipation of  excess  pore  pressures  in  the  vertical  direction 
as  well  as  in  the  horizontal  direction  parallel  to  the  cross 
section  of  the  embankment.   The  governing  equation  for  such 
a  condition  is: 


uihere 


c   =  coefficient  of  consolidation  in  the  vertical 

V 

d  irection 

c.  =  coefficient  of  consolidation  in  the  horizontal 

h 

d  irection 
u  =  excess  pore  pressure 
x<  z  =  the  horizontal  and  vertical  coordinate  directions 
t  =  time  after  the  excess  pore  pressures  were  created 

The  solution  for  equation  4. 18  depends  on  the  distribution 
of  the  excess  pore  pressure.   Therefore*  an  exact  solution 
is  not  possible  for  the  general  case.   An  approximate  solu- 
tion is  possible,  however*  by  specifying  the  distribution  of 
pore  pressure  on  a  grid  of  evenly  spaced  points  in  the  con- 
solidating layer  (Figure  4.4).   By  replacing  the  partial 
derivations  in  equation  4.18  with  finite  difference  approxi- 
mations on  this  grid,  it  is  possible  to  derive  the  following 
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FIGURE  4.4         FINITE  DIFFERENCE  GRID  FOR 

TWO-DIMENSIONAL  CONSOLIDATION 
EQUATION 


1A4 


finite  difference  solution  for  the  two-dinensional  consoli- 
dation equation: 

u.   .,^.=    a   Cu.^,   .,  +u.  ,   .,]       «C4.19> 

z    i.j+l,k     i,j-l,k 

+Cl-2a   -2a  Du.   ., 
X      z    1, J,k 

where 

i/j  are    column  and  row  identifiers  of  the  nodal  points 
on  the  grid  illustrated  in  Figure  4.4 
k  is  the  number  of  the  time  step 


c.  At  c   it 

h  V 


a      =  —  a      = 


^    <,ix>^  ^    <iz>^ 

ix.,ii!iz  sre    the  horizontal  and  vertical  spacings  of  the 
grid  shown  in  Figure  4.4. 

iit  is  the  time  increment  used  in  the  analysis. 

It  is  not  common  practice  to  evaluate  equation  4. 19 

because  c.      is    not  usually  measured.   Therefore.,  it  is  qen- 
h  - 

erally  assumed  that  consolidation  occurs  solely  due  to  vert- 
ical drainage.   The  governing  equation  for  this  condition  is 
the  uel 1-known  Terzaghi  one-dimensional  consolidation  equa- 
tion (Terzaghii  1943) 


£z 
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The  finite  difference  expression  for  the  solution  of  this 
equation  is 

u.  .   ,=a   u.  ,,+<!-  2a  >u.  ,+a   u.  ..<4.21) 
i,k+l     z   i+l,k  z    i,k     z   i-l,k 

The  accuracy  of  this  expression  is  ntaximized  uihen  Az 
approaches  zero  and  when  a   =1/6  < Perl off  and  Baron,  1976> . 

At  the  boundary  of  the  grid,  equation  4.21  must  be 
modified  to  account  for  drainage  conditions.   When  a  boun- 
dary is  drained*  the  excess  pore  pressure  is  assumed  to  have 
an  ambient  value  at  the  onset  of  consolidation  equal  to  half 
of  the  initial  excess  pore  pressure.   Thereafter!  the  excess 
pore  pressure  at  the  drained  boundary  is  set  equal  to  zero. 
To  calculate  the  excess  pore  pressure  at  an  undrained  boun- 
dary, it  is  necessary  to  assume  a  "mirror"  node  just  outside 
the  grid  with  a  value  of  excess  pore  pressure  equal  to  that 
of  a  node  just  inside  the  grid.   The  resulting  expressions 
for  nodes  at  the  top  and  bottom  boundary  respectively  are: 


u.  ,  ^,  =  2a   u.^,  ,  +  <1  -  2a  >    u.    .  <4.21a> 

1 ,k+l      zi+ljk  z    i,k 

and 


'-'•  ua.1  =  2a   u.  ,  .  +  CI  -  2a  )u.  ,  <4.21b> 

i,k+l      z   1-1, k  z    i,k 


When  the  consolidating  layer  is  composed  of  contiguous 
soil  layers  u»ith  different  values  of  c  ,    continuity  of  floui 
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must  be  satisfied  across  the  layer  interfaces.   Invoking 

Darcy's  law.  the  value  of  the  pore  pressure  at  the  layer 

interface  may  be  obtained  with  the  following  expression 
(Harr,  1966): 


i-H.k i-l.k  ,  .  22) 

■^i.k    ^i  +  l.k  ~  1  +  <k2/kj>/<A2j/iZ2>       '^•♦.^«i 


where 


k.^k   are  the  permeabilities  of  the  soil  above  and 
below  the  layer  interface,  respectively 
6z.,Az^  are   the  grid  spacing  above  and  below 
the  layer  interface,  respectively. 


As  a  first  approximation,  the  ratio  of  the  permeabilities 
between  the  upper  and  lower  layers  may  be  taken  to  be: 


where  c  .  and  c  _  are  the  coefficients  of  consolidation  in 
vl       v2 

the  upper  and  lower  soil  layers,  respectively. 


Once  the  excess  pore  pressure  is  calculated  at  all  the 
nodes  for  a  given  time  after  the  onset  of  consolidation,  the 
percent  consolidation  at  that  time  may  be  calculated  by 
evaluating  the  following  expression: 


H 

S   u.<z>dz 
z.       1 

UX  =  CI  -  ^ 1    X  100  <4.24> 

S   u<z)  dz 
0 
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inhere 

u.<z>    =   the    initial    distribution   of   excess   pore 

1 

pressure  with  depth 
u<z)  =  the  distribution  of  excess  pore  pressure  with 

depth  at  the  time  in  question 
H  =  the  thickness  of  the  consolidating  layer 
z  =  coordinates  in  the  direction  of  the  depth 

Equation  4.  24  can  be  used  to  calculate  the  percent  con- 
solidation in  each  layer  of  a  series  of  continguous  soil 
layers  as  well  as  the  overall  percent  consolidation.   A  pro- 
gram for  performing  these  computations  is  provided  in  Appen- 
dix H. 

Values  of  the  initial  excess  pore  pressure  that  are    to 
be  used  in  the  analysis  should  be  calculated  with  equation 
4.9.   This  insures  that  the  distribution  of  excess  pore 
pressures  used  to  calculate  the  percent  consolidation  is  the 
same  as  the  distribution  used  to  calculate  the  magnitude  of 
consolidation  settlement. 

The  coefficient  of  consolidation  is  defined  by  the  fol- 
lowing relationship: 


c   =  r <  4  .  25  ? 

V     a   0 
V   w 


where  y   is  the  unit  weight  of  water, 
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Since  the  values  of  the  permeability*  the  void  ratioi 
and  the  vertical  compressibility  of  the  soil  change  as  the 
consolidation  progresses*  it  is  necessary  to  simplify  the 
characterization  of  the  value  of  c   during  consolidation. 
This  can  be  done  with  the  controlled  gradient  consolidation 
test  which  is  used  in  developing  the  preconsol idation  pres- 
sure profile.   The  value  of  c   obtained  in  a  controlled  gra- 
dient test  is  calculated  with  the  follo-t;ing  expression 
(Lowe/  Jonas  and  Obrician.  1969): 


c=7rr—  <4.2fe^ 

V    &t    2u 


where 
6c 


=  time  rate  of  change  of  applied  stress 
H  =  the  sample  thickness  ' 

u  =  excess  pore  pressure  maintained  at  the  undrained  end 
of  the  sample 

Thiis  expression  allows  the  calculation  of  c   continuously 
during  consolidation  without  recourse  to  either  the  loga- 
rithm or  square  root  of  time  curve  fitting  methods.   Other 
advantages  of  this  testing  procedure  include: 

1>  the  test  may  be  run  at  low  strain  rates  that  ap- 
proach consolidation  rates  in  thr  field 

2>  the  excess  pore  pressure  is  appro>:imately  constant 
across  the  sample 

3>  secondary  compression  does  not  oicur 
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Once   the   value    of    c       of    a   soil    layer    is   Jcnouir*   along   a 

V 

compression  curve  like  the  one  shouin  in  Figure  4.  la>  it  is 
possible  to  choose  a  representative  value  of  c   correspond- 
ing to  the  average  value  of  vertical  effective  stress  during 
consolidation.   As  was  the  case  for  calculating  the  magni- 
tude of  settlement,  this  procedure  should  be  repeated  in  a 
number  of  artificial  strata  within  the  soil  layer  because 
the  soil  parameters  and  the  stress  changes  will  vary  with 
depth. 

Budget  and  time  constraints  may  prohibit  the  type  and 
amount  of  testing  necessary  to  perform  the  analysis  just 
described.   When  this  is  the  case,  an  estimate  of  the  varia- 
tion of  c   tuith  depth  needs  to  be  made.   In  general.,  the 
value  of  c   in  a  soil  that  is  overconsol i dated  will  be  sub- 

V 

stantially  higher  than  c   of  a  normally  consolidated  sample 
of  the  same  soil.    If  it  is  possible  to  estimate  which  por- 
tions of  the  soil  layer  will  be  normally  consolidated  and 
overconsol idated  during  the  consolidation  process, 
the  entire  layer  can  be  reduced  to  a  two-strata  sys- 
tem.  Normally,  c   of  the  upper  portion  will  correspond  to 

an  overconsol idated  state  and  c   of  the  lower  portion  will 

V  ^ 

correspond  to  a  normally  consolidated  state. 

Example  4.  3 

It  is  desired  to  estimate  the  time  rate  of  settlement 
of  the  consolidation  settlement  at  the  centerline  of  the 
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embankment  in  Example  4.1.   The  liquid  linit  of  the  soil  in 
the  coTisol  idat  iT»9  layer  is  45.   Typical  values  of  c   are 
8.61  m^'/day  CO.  8  ft^/day)  and  2.15  m^/daj  <0.2  ft^/day>  for 
the  overconsol idated  and  normally  consolidated  portions  of 
the  layeri  respectively  (Holtz  and  Kovacsi  1981).   The  ini- 
tial excess  pore  pressures  mere  estimated  using  equation 
4.  10.   The  values  of  the  stress  changes  uere  taken  from 
Table  4. 1.   It  was  assumed  that  the  overconsol idated  portion 
occupies  the  upper  6.  Im  (20  ft)  of  the  soil  layer.    In  this 
example  the  bottom  of  the  consolidating  layer  is  undrained. 
The  results/  obtained  with  the  computer  program  in  Appendix 
H/  are    shoun  in  Table  4.5.   The  time  rate  of  settlement  in 
the  normally  consolidated  stratum  proceeds  at  a  much  slou/er 
rate  than  the  rate  of  the  entire  layer.   Since  most  of  the 
settlement  occurs  in  the  normally  consolidated  stratum,  it 
would  be  unconservat i ve  to  estimate  the  time  rate  of  settle- 
ment on  the  basis  of  the  overall  rate  of  consolidation  for 
the  soil  layer.   A  better  estimate  of  the  time  rate  in  this 
case  can  be  made  by  assuming  that  all  of  the  settlement 
occurs  in  the  normally  consolidated  stratum  and  proceeds  at 
the  rate  computed  for  this  stratum.   This  calculation  is 
included  in  Table  4.5. 
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V  -  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  WORK 

Cone lusions 

1.  A  hybrid  method  of  compaction  specification  which 
makes  the  insitu  water  content  of  the  embankment 
soil  equal  to  the  OMC  is  introduced. 

2.  The  Simplified  Bishop  factor  of  safety  option  in 

the  STABL  program  has  been  recoded  to  correct  various 
difficulties. 

3.  A  methodology  for  adjusting  the  Simplified  Janbu 
factor  of  safety  that  is  used  by  several  STABL 
options  to  definitions  of  the  factor  of  safety  that 
are  more  familiar  is  presented. 

4.  Embankment  side  slope  design  has  been  illustrated 
for  short  and  long  term  situations  using  laboratory 
compacted  shear  strength  data. 

5.  Geometric  and  probabilistic  interpretations  of  the 
factor  of  safety  are    introduced  to  illustrate  their 
usefulness  in  the  selection  of  an  appropriate  factor 
of  safety  for  design. 

6.  A  methodology  of  predicting  settlement  of  embankment 
foundations  has  been  illustrated.   Co-nputer  programs 
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to  compute  the  magnitude  and  time— rate  of  settlement 
are    included.   These  programs^  in  conjunction  with 
STABLj  form  an  analysis  package  for  the  design 
of  embankments. 

R  e  c  ommen  d  a  t  i  on  s 

1.  To  use  the  hybrid  approach  for  specifying  compac- 
tion>  values  of  the  coefficient  of  compaction 
that  reflect  the  influence  of  soil  type*  compactor 
type>  the  operating  procedure^  and  the  number  of  passes 
need  to  be  developed. 

2.  The  current  version  of  STABL  should  be  augmented  by 
adding  a  complete  equilibrium  method  of  calculating 
the  factor  of  safety  such  as  the  Spencer  method 
(Spencer,  1973). 

3.  The  LEMIX  program,  uthich  calculates  the  factor  of 
safety  on  a  three-dimensional  surface,  demonstrates 
the  importance  of  three-dimensional  effects.   The 
next  step  in  this  field  should  be  the  development 
of  a  program  that  allous  the  random  generation  of 
general  shaped  three-dimensional  surfaces  as  uell  as 
the  subsequent  calculation  of  the  factor  of  safety  on 
these  surfaces. 

4.  The  geometric  and  probabilistic  approaches  uhich  were 
introduced  should  be  further  developed  to  determine 
acceptable  design  values  of  the  geometric  factor  of 
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safety  and  the  probability  of  failure  for  use  in  slope 
design. 
5.   Procedures  should  be  developed  to  apply  the  statisti- 
cal relationships  between  the  material  parameters  of 
lab  and  field  compacted  clays  in  settlement  predictions 
within  an  embankment  and  in  slope  stability  analysis. 
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APPENDIX  A 

DERIVATION  OF 
THE  SIMPLIFIED  BISHOP  FACTOR  OF  SAFETY 


Notation  after  Siegel  (1975). 

STEP  1  -  Enforce  moment  equilibrium  of  sliding  circular 
mass  divided  into  n  slices: 


Sh   =  0 
o 

ii 

1  ^                ^ 

n  T\ 

-    S    CiS    :  R    -    S  C<AU^    sinje    +    AQ    siT»6>CRcoso:    -    h)D 

1            ^  1               ^ 

H    E    CK     6W<Rcosoc    -    h       >D  <A.1> 

J         h  eq 

u>here 

R  =  radius  of  the  circle 

i.S   =  ~  CC-   +  iN'tan*--  1  <A.2> 

r    F5     a  a 

C''   =  c  ■  dx/cosa 


Dividing  equation  A.  1  by  R  yields: 
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o  .  V  £ 

Ti  ri 

-    2  [.iS    D    -    E    t  <£.LI     '  sin.e    +    AQ    sirii;":;coso:    -    h/R> 

r.  h 

+    E  Ck.     .iW -Ccoso: 1^>  D 

J  r  I  K 


STEP  II  -  Substitute  iS   into  equatiori  h.3  and 

assume  the  factor  of  safety  is  equivalent  in  each 
slice: 


n 
S 

1 

<C'       + 

a 

i>N  ■■  t  an*  ■'     > 

a 

n 
Z 

1 

n 

%   *   f   ^5 

FS  = : '::A.4> 


iL'here 


A.-,    =    t: liU    I.' 1 -k     >    +    liU  ^  ■  cos/J    +    liQ  '  co=-t>    sino:  ■:  A    4a  > 

A  .    =    t: liLI ^,    s i nj£!    +    i.Q    si n .£ >  •:: c c^o:    -    h/k ':■'  <'.  A  .  4bi  :■ 


A^    =    k,     liNtcosa    -    — ^.l  CA    4r  1 


STEP    III    -    Sum    forces    in    \yertical    direction    for    each    slice 


iWCl-k^^    -    CC-'.    +    £,N-'tan+'    :i'sincc/FS    -    iN'cosc- 
V  d  a 

+    AQ'cosS    +    ALkcosje    -    iU      coscc  CA    F, 


Rearranging    equation    A. 5   yields 


AN  ■'■':;  t  an+ •■'       sincc/FS   +    coso.'>    = 

a 

^Wt:i-k     >    -    C ■•■  _  ■  s i na/FS    +    AQcos.S    +    ■iU,cos.t:    -    iU      coscr    (A    A> 
V  •=!  .t!  a.  ■  " 
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Final  ly. 


AWCl-k    y-C"      siTia/FS+<iQcos«+6U^coi:je-iiU      cosa 
~  cosa    +   tan*"      sina/FS  Kf< .  ,.> 


Substituting    equation    A. 7    into    equation    A. 4    yields: 


FS    =    CC       +    tan*--     ■<iW<l-k    >    -    C       sincc/FS 
a  a  V  a 

+    iQ    cos5    +    AU^cos^    -    AU    ■cosa>/ 
JO  a 

Ccoso:    +    tan*"     '  sinoc/FSD/CSA^    -   SA.    +   SA_>  <A.8> 

a  ^'4  0 


or  rearranging 


n    C    +tan+''    seco:<AWCl-k     >+AQcos£+AU,cos;e-6U   cosa> 
2[;_-§ § V fi oc 

.  l+tan*-'       tana/FS 

p-_         1 a 

FS    = — <A.9> 

Ti  n  n 

E    A-    -   S    A .    +    E    A_ 

1      ^         1       '^         1       ^ 


For    simplicity    of    coding    equation    A.  9    may    be    written    as 
f ol lows : 


n                 A. 
S    [■ ~ 


J       1    +    A^/FS 

PS    = CA.IO:- 

n  n  n 

S   A-    -    S    A.    +   S   A^ 

1      ^         1      '^         1       ^ 


where 


•^1    =    C:'_    +    tan-*-'       secoiCAWCl-k     >    +    AQcos6  <:A.10a> 

X  44  3  V 


+   ALl.cosje   -    AU      cosoc? 
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A   =  tan  *•'   tana  CA.10b> 

This  expression  is  programmed  in  STABL3. 

The  changes  in  the  STABL  code  necessary  to  program  equation 
A  10  are    included  in  the  follou/ing  pages. 
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APPENDIX  B 
CALCULATION  OF  THE  RADIUS  OF  A  SPECIFIED  CIRCULAR  SURFACE 

STAEL  randomly  generates  circular  surfaces  by  generat- 
ing successive  chords  of  a  circle  which  are    inclined  at  a 
deflection  angle  with  one  another.   Since  all  of  the  chords 
are    of  equal  length*  the  radius  may  be  determined  with  the 
follouiing  equation: 

R  =  5  sir,  <^>  a;.l> 

where 

R  =  radius  of  the  circle 

T  =  chord  length  of  the  segments  circumscribing  the  circle 
lJh  -  deflection  anqle  of  the  segments  c  i  rcLnTi;  c  t  ibi  tiq  tf-.e 
circle 

However,  when  a  specified  surface  is  input,  the  chord 
lengths  and  deflection  angles  will  not  have  constant  values. 
In  this  case  a  more  general  expression  for  calculating  the 
radius  is  required.   To  do  this  any  two  adjacent  segments  on 
a  known  circular  surface  may  be  taken  (Figure  E. 1). 
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The    equation    of    line    PI    which    is    the    perpendicular 
bise-ctor    of    the    chord    defined    by    '^x    ,y..j    and    <v.., ..  y..,  >    is 


u  =  -  z r  "•'  x-x   .  '■   +  u   .  <  B  .  2 . 

^.y^-y^.'  cl  -cl 

u;here 


<X.  ^y^>  are  t^1e  coordinates  of  the  first  point  or  the 

first  chord 
Cx,,,  y..,;j  are  the  coordinates  of  the  second  point  on  the 

first  chord  and  the  first  point  on  the  second 

chord 

'.x  .  ..  u   ■;•  ^m    the  coordinates  of  the  intersection  of 
c  1  "  c  1 

line  PI  and  the  first  chord. 

Similarly,  the  equation  of  line  P2,  which  if,  the  perpendicu- 
lar bisector  of  the  chord  defined  bu  <x..  ,u.-,  •'  and  ';x..  ,u.->..  i; 


u  =  —  (.  X— X  ^  :>    +  u  ,  (■  B 

="      yo-y-:,      c2     ='c2 


tlx.-,,  y-, :■  are  the  coordinates  of  the  second  point  on 
the  second  chord 
'^X  ... ..  u  ,/'  are  the  coordinates  of  the  intersection  of 

C  CL   "  C  c! 

line  P2  and  the  second  chord. 

Lines  PI  and  P2  intersect  at  the  center  of  the  circle. 
The  coordinates  of  the  center  of  the  circle  may  be  found  by 
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setting  equations  B.  2  and  B.  3  equal  to  each  other.   The 
resulting  expression  is: 


^-^<x  -- V^)  +  -hr^    =  -2— ^<x  -  -^^>  +  -^r-^  <B.4> 


y^-yj   o    2        2       yg-yg   o     2 
Rearranging^ 

x^  -  c<Xj-x2>  -  <x2-x^><y2-yj>  +  <y3-gi><y2-yi><y3-y2>3/2/ 

tCXj-XgXy^-yg?  -  Cx^-XgXyg-yj)^  <B.5> 


The  value  of  y   may  be  obtained  by  substituting  the  value  of 
X   in  equation  B.3.   The  result  is 

X  ^x 
y   =  -^ — ^  <x  -X  ->  +  y  ^  <B.6> 

^o   yg-ya   ©  c2    «c2 

where 


^c2  =  <>«2***3^^2 

«c2  -  ^y2-^y3^/2 


Finally,  the  radius  may  be  found  with  the  following  expres- 
sion: 


^   =  N^'*o-^2>^  *  ^yo-«2>^  ^^-^^ 


These  equations  Br&   coded  in  the  current  version  of  the 
STABL  program . 
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APPENDIX  C 
THE  FRICTION  CIRCLE  FACTOR  OF  SAFETY 

The  stability  number,  c/<FSVH>,  of  a  trial  circular 
surface  through  a  homogeneous  slope  such  as  the  one  shoun  in 
Figure  C.  1  may  be  determined  by  the  following  sequence  of 
calculations  <Tay lor> 1940) : 

n  =  ^  <cotx  -  coty  -  cotB  +  sin*' cscx 'cccy)     <C.1> 
where 

&   =  sides lope 

♦  =  friction  angle  of  the  slope 
x,\i    <=   angles  shown  in  Figure  C.  1 

If  n  >  0<  the  trial  circle  passes  under  the  toe  of  the  slope 
as  shown  by  surface  III  in  Figure  C.  2.   If  n  <  O.  the  trial 
surface  starts  at  the  toe  as  shown  by  surface  I  or  II  in 
Figure  C.  2. 

When  n  <  0/  equations  C. 2  through  C. 5  are    evaluated. 

12        2 
J,      2  ^^^  ***^y  ^^^   y— coty>  +  cotx  -  cotB 

—  = _ <C.2> 

2<l-2cot    j8>    +   cotje<cotx-coty>    +   cotx'coty 
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H  2 

coto  «   r4y 'Eecx'cscx' CSC    y   -   tanx  <C.3> 


siTi<u-v>    =   r-rsin   ucscx   cscy 'sin*  <C.4> 

eiO 


12  2 

-  esc    xCycsc    y   -   coty>    +   cotx-cotB 
■^rr:-  =   ^ : — : : r <C.5> 


FS   V   H  2cotxcotv    +   2 

When    n   >   O.     equations    C.  6   and    C.  7   replace    equations    C.  2   and 
C.  5. 


zg  =    c|  csc^x    <ycsc'^y-coty>    +   cotx   -   cotiB   -   2t»3/ 

1  2 

C^<l-2cot    B>    +   cotB' <cotx-coty>    +   cotxcoty 

+   2n      -   2n   sin* 'cscxcscyD  <C.6> 


12        2 

—  CSC  x<ycsc  y  -  coty)  +  cotx  -  cot*  -  2n 

— — <C.7> 


FS  y  H  2cotx  cotv  +  2 

Given  B   and  ^,     it  is  possible  to  search  for  the  surface 
on  uihich  the  stability  number  has  a  maximum  value  bg  check- 
ing all  combinations  of  the  angle  x  from  0  up  to  B  and  all 
the  values  of  the  angle  y  from  0  to  90  .   The  search  may  be 
limited  to  y  <  x  if  the  limit  equilibrium  surface  does  not 
go  beneath  the  elevation  of  the  toe.   The  value  of  the  sta- 
bility number  so  obtained  is  almost  identical  to  that  which 
is  presented  in  the  stability  chart  that  Taylor  developed 
for  simple  slopes  in  homogeneous  material  except  that  the 
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equations  presented  here  do  not  reflect  the  correction  for 
the  error  associated  uiith  the  assumption  that  the  resultant 
of  the  normal  and  frictional  forces  acts  tangent  to  the 
friction  circle.   The  maximum  value  of  this  correction  is 
approximately  7%  (Taglor<  1937). 

To  calculate  the  factor  of  safety  of  a  slope,  the  fac- 
tor of  safety  must  be  the  same  on  the  cohesion  intercept  and 
the  friction  angle,  i.  e.  . 

FS   =  FS.  «CC.8> 

c      •♦ 

This  condition  may  be  satisfied  uith  the  following 
iterative  procedure  uihich  is  used  by  the  computer  program 
included  hereafter: 

1 .  Assume  a  value  of  FS^ 

2.  Calculate  the  required  value  of  ♦  as  follows: 

♦     =  tan"^  <^f5^>  (:C.9> 

req  FS^ 

3.  Calculate  the  maximum  stability  number^  N  ,    for  ♦ 

*  s       req 

with  equation  C.7. 

4.  Calculate  FS   mith  the  expression: 

FS   =  <:c/VH>/N  <C.10> 

c  s 

5 .  Compare  FS   and  FS . 

a>  If  *  =  0,  then  FS  =  FS  . 

b>  If  ♦  9*  0  and  FS.  =  FS  ,  then  convergence  has 

been  obtained 
c>  If  ♦  ?«  0  and  FS^  ^   FS^,  repeat  steps  2 

through  5. 
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This  procedure  luorks  except  when  the  stability  numbers  close 
to  zero  cause  numerical  instability. 
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User  Manual  -  Friction  Circle  Factor  of  Safety  Program 

The  following  program  uias  developed  in  77  Fortran  on  a 
CDC  6000  series  computer.   All  input  is  in  English  units. 
Houever,  any  d imensional ly  homogeneous  set  of  units  will 
work.   All  input  is  unformatted. 

1.  Input  on  one  record: 

a)  the  friction  angle  of  the  slope  (degrees) 

b>  the  slope  angle  (degrees) 

c)  the  cohesion  intercept  (psf) 

d)  the  slope  height  (ft) 

e)  the  density  of  the  soil  (pcf) 

2.  Input  on  one  record: 

a)  an  integer  variable  controlling  the  limits  of 
the  search: 

if  the  slip  circle  is  completely  above  the  elevation 
of  the  toe  input  '1' 

if  the  slip  surface  intersects  the  toe  of  the  slope  and 
may  descend  down  to  a  specified  depth  limit<  input  '2' 
if  the  slip  surface  may  pass  under  the  toe  of  the 
slopei  input  '3' 

3    If  the  value  of  the  integer  variable  on  record  #2  is  '2'i 
input  on  one  record: 
a)  the  depth  factor  of  the  search. 

The  depth  limit  is  the  depth  beneath  the  top  of  the 
slope  to  the  bottom  of  the  deepest  circle  that 
is  geometrically  possible. 
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The  depth  factor  is  obtained  by  dividing  the 
depth  limit  by  the  slope  height.   The  depth  limit. 
DH.  is  illustrated  in  Figure  D.  2. 

If  the  value  of  the  integer  on  record  #2  is  '3'. 
input  on  one  record: 

a)  the  depth  factor  of  the  search. 

b)  the  toe  factor  of  the  search. 

The  toe  factor  is  the  multiple  of  the  slopes  height 
that  the  slip  circle  can  extend  beyond  the  slopes  toe 
(  see  Figure  D.  2) . 


The  program  listing  is  included  in  the  following  pages. 


ISA 


c  a  mm  on  /  hate/  num.X;  numy;  toe.  slope,  dangle. 
1  d ep th f , n 1 imi t , nsmax 

;Dmmon/lo5t/  coti/ sinphi 
re-al  n5max,nli(nit 
integer  toe- cycles 
data  toll,  tol2*  cycles,  fsphi/l.  Oe-8,  G.  01,  100,  1.  0/ 


pi  =  a  c  0  5 ( - 1 .  0 ) 


10 


12 


14 

16 

c 
c 


read  (5,  *)  phi,  slope,  cohes-hiqht,  gamma 
lurite  (6,10)  phi,  cohes,  gamma,  hight,  slope 
■Format  (  '  friction  angle=',  f6  1,  'degrees',  /, 


1  '  cohe5ion=  ',  fS.  1 ,  /, 

2  '  density=',  f6.  1,  /, 

3  '  hight=',  f6.  1.  /, 

4  '  side-sl  ope= ',  f  6.  1,  ' 
read  (  5,  *  )  toe 

if  (toe.eq.  1)  u)rite<6,  12) 
format  (  '  circles  start  at  the 


if  ( toe.  eq.  2) 
if  ( toe.  eq.  3) 
i  f  ( toe.  eq.  2) 

format  ( '  D  = 
if  (  toe.  eq.  3  ) 
format     ( '     D    = 


veati  (  5,  «)    depthf 

read(5,*)    dep  th  f ,  nl  i:-ni  t 
ujrite<6,  14)    depthf 

',  f5.   1,  //) 

write(6, 16)  depthf, nlimit 
',  depthf,/,  '  n  =',nlimit,/', 


degrees ' ,  / > 


toe  w'lth  depth  factor  =  l 


//  ) 


numx  =  intCslope)  -  1 

numy  =  39 

if  (toe  eq. 1)  numy  =  numx 

slope  =  slope*p 1/180.  0 

pn 1  =  ph i*p  i/iaO.  0 

if  (abs ( 5lape-pi/2.  0)  1 1.  tol 1 )  then 

coti  =  0.  0 
else 

coti  =  1. 0/tan ( slope ) 
end  1  f 

dang le  =  pi/180. 0 
sinphi=sin(phi) 
cgamh  =  c  ohes/ (  gamma-K-h  i  g  h  t ) 


do  100  1  =  l^cycles 

phireq  =  atan < tan ( ph i ) /f sp h i ) 

sinphi  =  sin(phireq) 

call  p  h  1  c  1  r  . 

fscohs  =  c  gamh/nsma.x 

if  (phi.lt.  tol2)  then 

print*,  'fE=',  fscohs 

print*,  'stability  numb er= ' ,  nsma x 
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ioo 

c 
c 

c 
c 
c 


stop 
e  n  d  1  f 

fsdiff  =  abs<f5coh5-f5phi) 
pr  int*)  f  sp  h  i/  f  5C  oh  5.  nsma  X 
:f  ( fsdiff . It.  tol2)  then 

print*-  '  f5=  '.fscohs 

stop 
elseif  ( fsdiff .  ge.  tol2)  then 

fsphi  =  <f5phi  +  fscohs)/2.0 
end  if 
if  (i.  eq  cycles)  then 

print*,  '  convergence  not  obtained 

stop 
end  i  f 
ccntinue 


end 


subroutine  phicir 

CO mm on/hate/  numx*  nurny.  toei  slopes  dangle, 
1  dep th f > nl imi t< nsmax 

common/lost/  coti,5inphi 

c  OfTimon/ 1  ove/  en,  x,  y,  cscx,  cscy,  cotx,  cotu^  stbnum 
■i^eai  nsmaxinlimit 
int eg  er  toe 
y  =  0.  0 
nsma  x  =  0.  0 


do 

y  = 

cot 

CSC 
X  = 
do 

X   = 

If 

cot 

CSC 

en 

i  ■? 
i  f 
d  = 

i-  J  1 


C 

c 


25 

y 
y  = 
y  = 

0. 
25 

X 

(to 
X  = 
X   = 

=  0 
(  <  t 
f  to 
0. 
(to 

1   5 

a  X 
tin 


er y  =  1 , numy 

dangle 

.  0/tan<  y ) 

.  0/sin ( y  ) 


er  X  = 
dang  1 
eq.  1. 
.  0/ta 
.  0/si 
*(cot 
.  eq.  1 
eq.  3. 
(cscx 
ne.  1 . 
bnm 
ama  x  1 


x.lt.y)     go    to    25 


1  <  ni 
e 

and. 
n(  X  ) 
n(  X  ) 
X  - 
.  or. 
and . 
•K-csci 
and.  d.  gt.  depthf.  and.  y    gt.  x>    go    to    25 


cotu  -  coti  + 
toe.  eq.  2  ) .  and 
en.  g  t.  nl  imi  t ) 
:y  -  cotx*coty 


sinphi*cscx*cscy ) 
en.gt.O.  0)  go  to  : 
go  to  25 
+  1.  0) 


( nsma  x ,  stbnum ) 
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end 
c 
c 
c 

subroutine  stabnm 

common/love/  en<  x, y, cscx,  cscu- cotx/ coty, stbnum 

CO mm on /lost/  cotijsinphi 
c 
c 

c      determine  the  stability  number  of  a  specified  slip 
c      surface  on  a   specified  slope. 
c 
c 
c 

sec  X    =  1.0/coB<x) 

1  f  (  en.  1  e.  0.  O)  then 
c      toe  circle 

paraml  =  (0.  5-«-c  sc  x»*2*(  y*c  sc  y**2  -  coty)  +  cotx 

1  -  coti  )/<  1.  0/3.  0*<  1.  O  -2.  0*coti**2)  +  coti 

2  «(cotx  -  coty)  +  cotx*coty) 
param2  =  paraml*y  *sec  x*c  sc  x*c  sc  y  ■»*-<-2  -  tan  (  x  ) 
u  =   atan ( 1 . 0/param2) 

param3  =  paraml*sin < u ) *c sc x*c scy *s inph i 
uv  =  asin(param3) 

V  =  u  —  uv 

Btbnum  =  (  0  5-fi-c  sc  x-«-«-2«-(  y -^-c  sc  yK-*2  -  coty^  +  cotx  - 
1  coti  )/ <2.  Oifcotx/tan(v)  +  2.0) 

elseif  (en.  gt.  0.  0)  then 
c      slip  surface  beneath  the  toe  of  the  slope 

paraml  =  (  0.  5*c  sc  x-M"K-2*(  y  *c  sc  y  ■«-*2  -  coty)  +  cotx 

1  coti  -  2.  0*en)/(l.  0/3.  0*(  1.  0  -  2.  0*coti 

2  **2)  +  coti»<cotx  -  coty)  +  cotx-R-coty  + 

3  2.  0-«-en**2  -  2.  0-»-en-»-s  inp  h  i*c  sc  x*c  sc  y  ) 
param2  =  par  ami  *y  #sec  x -o-c  s  c  x -fi-c  sc  y  *-*2  -  tan(x) 

u  =  atan(l.  0/param2) 

paramS  =  paraml  *s  in  (  u  ) -s-c  sc  x  *c  sc  y  «s  i  np  h  i 

uv  =  asin(param3) 

V  =  u  -  uv 

stbnum  =  (  0.  5-«-c  sc  x-»-#2»  (  y«-c  scy  »*2  -  coty)  +  cotx 
1  coti  -  2.  0*en  )  /  (2.  0»cotx/tan  (  v)  +  2.0) 

end  1  f 
return 
end 
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APPENDIX  D 
THE  BETA  DISTRIBUTION 


The  beta  distribution  for  the  variable  x  is  defined 
by  the  following  density  function  (Harr.  1977): 

^ ,    .         _1 r    Ca+j6+2>  ,x-a.  PC    .b-x.jB  ._    .  . 

f<x>    =  r —  :pr~: — . . .    -    .  .    , .    <r — >      (.r — ^  <D.1> 

b-a   r    <a+l>    r    <;e+l>      b-a  b-a 

inhere 

b  =  upper  bound  of  the  density  function 
a  =  louer  bound  of  the  density  function 
r  =  Gamma  function 

The  expected  value  of  the  variable  x  is 

and  the  variance  of  the  variable  x  is 

y<x>  =   <b-a>^  <«  I''    '^*    '>  <D.3> 

"Ca  +  jS  +  2>   Ca  +  jS  +  3> 

Uhen  a,     b>  E(x}  and  V(x)  are  known,  the  constants^ 
a  and  £,    may  be  defined  by  the  following 


^2 
a   =   ^^   <1    -   Si:>    -    <1    ■*■   Si>  <D.4> 


B   =  SLJt_l  _  <(x  +  2>  <D.5> 


where 


a  = 


E<x>    -   a 
b   -   a 


^   =         V<?<> 


<b   -   a>' 
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<D.6:» 


<D.7> 
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APPENDIX  E 

PROBABILISTIC  SLOPE  STABILITY  AMALYSIS  BY 
THE  POINT-ESTIMATES  METHOD 


The  folloujing  program  was    developed  to  determine  the 
probability  of  failure  of  a  simple  slope.   Failure  is 
defined  to  be  a  condition  uhen  there  is  any  surface  along 
uhich  equilibrium  cannot  be  maintained.   Limit  equilibrium 
is  checked  with  the  Taylor  friction-circle  method.   The 
values  of  the  soils  friction  arigle,  ^,    cohesion  intercept., 
c,  and  density.,  "^ ,    are  assumed  to  be  syir.Tictr ical ly  distri- 
buted.  The  mean  value  and  the  standard  deviation  of  the 
strength  factor  of  safety  of  the  slope  aie  obtained  with  a 
tuo-point  estimate  (Rosenb 1 uethi 1975) .   A  beta  distribution 
is  fitted  to  these  statistical  moments  and  the  upper  and 
louer  bounds,  of  the  factor  of  safety.   The  probability  of 
failure  is  computed  by  numerically  integrating  equation 
3.  13. 
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User  Manual  -  Probabilistic  Slope  Stability  Program 

This  program  luas  developed  on  a  Vax  computer.   All 
units  are    in  English  units*  but  any  d imens i ona 1 ly  homogene- 
ous set  of  units  will  work.   All  input  is  unformatted. 

Parameters  Defining  the  Distribution  of  the  Values  of  the 
Material  Properties  of  the  Slope 

1.  Input  on  one  record: 

.=!>  fTi'E'aTi  value  of'  4-  <C degrees > 

b>  the  coefficient  of  variation  of  •+ 

c>  the  louier  bound  of  4-  iC degrees > 

d>  the  upper  bound  of  4-  t! degrees > 

2.  Input  on  one  record: 

a)  the  mean  value  of  c  <psf) 

b)  the  coefficient  of  variation  of  c 

c)  the  loujer  bound  of  c  (psf) 

d)  the  upper  bound  of  c  (psf) 

3.  Input  on  one  record: 

a.:'  the  mean  value  of  V  <pcf > 

b>  the  coefficient  of  variation  of  V 

c  ■.:■  the  lomer  bound  of  V  <pcf> 

d;:-  the  upper  bound  of  y  <pcf> 

Parameters  Defining  the  Geometry  of  the  Slope 

4.  Input  on  one  record: 
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a)  the  height  of  the  slope  (ft) 

b)  the  slope  angle  (degrees) 

Parameters  Defining  the  Type  of  Limit  Equilibrium  Surface 

5.   Follow  instructions  2-4  in  the  User  Manual  for  the 

Taylor  friction  circle  factor  of  safety  which  is  found  in 
Append  i  X  C. 


The  program  listing  is  included  in  the  following  pages. 
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5 
8 

12 

14 

16 

c 
c 
c 

17 


c ommon /gr oovy /  n 1 i mi t» d ep t h f 

common/hate/  numx<  numy,  toei  dangle/  nsmax 

real  mean(3)*  5tddev(3)>  cv<3)/  bound(3/ 2))  nlimit, nsmax 

integer  toe 

data  f  5x,  iout  /I.  0,   1/ 


input  material  parameters 
a)   mean  value     b)   coe-Pficient  of  variation 
c)   louer  and  upper  bounds  of 

1)  the  friction  angle 

2)  the  cohesion  intercept 

3)  the  soil  density 


do  5  i  =  1.  3 

read  (5i*)  mean  (  i  )»  cv(  i  )/  bound  <  ii  1  ).  bound  (  i  *  2) 

stddev(i)  =  mean ( i )«cv< i ) 

continue 

ujr  1  te  (6.  8) 

format  (15(/)/t46»  'soil  propert i es  ' ,  /// 

1  t20,  'mean',t40i  'coeff.  of  var.  ' ,  thO,   'lower  bound 'i 

2  t80<  'upper  bound '<//) 

write  (6,12)  mean ( 1 ), cv < 1 ), bound ( 1 , 1) , b ound < 1 , 2 ) 
format  ('  p  h  i  ' ,  tlS>  f  6.  2,  t42..  f  5.  3,  t62<  f  6.  2,  tB2,  f  6.  2  ) 
urite  (6,14)  mean (2) , cv ( 2 ) , b ound ( 2, 1 ) , b ound ( 2, 2 ) 
format  ('  cohesion',  tl6,  fS.  1,  t42,  f5  3,  t60,  f8.  l,tSO,  fS.  1) 
write  (6,16)  mean ( 3) , c v ( 3 ) ,  b ound (3,  1 ) ,  b ound ( 3, 2 ) 
format  ('  densi  ty ',  1 19,  f  5.  1 ,  t42,  f  5.  3,  t63,  f  5.  1 ,  t83,  f  5.  1,  5(  /  ) , 
1         t45,  '  slope  geometry',//) 


input  geometry  of  the  slope 


=  ',  f5.  1,  /,  t37,  '  slope  =',  f5. 


read  (5,*)  hight,  slope 
write  (6,17)  hight, slope 
format  (t37, '  slope  hight 
numx  =  int(slope)  -  1 
numy  =  89 

pilSO  =  acos(-l.  0)/180.  0 
dangle  =  pilSO 
slope  =  s lope*p i 180 


choose  the  type  of  failure  surface 

if  the  depth  factor  =1.0 
if  the  depth  factor  >  1.0 

if  the  failure  surface  passes  under  the  toe 
if  toe  option  =  2     input  the  depth  factor. 
ir    toe  option  =  3     input  the  depth  factor 

nl imi  t 


set 

toe  =  1 

set 

toe  =  2 

set 

toe  =  3 

dep thf 

and  the 

toe  factor 
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'ead  (  5,  ■»■ )  toe 


18 


19 


c 
c 
c 

21 


23 


25 


if  (toe.eq.  1)  then 

ur i te  (6. IB) 

format  (t37;  '  circles  start  have  depth  f ac tor=l .  0  ' .  // ) 

numy  =  numx 

depthf  =  1.0 

nlimit  =  1.0 
end  i  f 
1  f  ( toe,  eq.  2)  then 

read  (5, *)  depthf 

ujr  1  te  (6.  17)  depthf 

format  (t37.  '  circles  start  at  the  toe  luith  depth  factor=',f6. 
1  //) 

nlimit  =  1.0 
end  if 
if  ( toe.  eq.  3)  then 

read(5i*)  dep  thf  >  nl  imi  t 

write  (6>20)  depthf>nlimit 

format  (t37j  '  circles  extend  beneath  the  toe'//. 

1  '  depth  factor  =',f6.  1./, 

2  '  toe  factor    ='if6.  1.//) 
end  if 

calculate  the  value  of  the  central  factor  of  safety 

'ojrite  (6,21) 

format  (  1  h  1.  10( / ) ,  t20,  '  ca 1 culat i on  of  central  f s  '  ) 

phmean  =  mean  ( 1  ) -K-p  i  ISO 

call  fsafty  ( s 1  ope,  h i gh t,  phmean< mean (2 ) < mean ( 3 ) ,  f scntr  ) 

lur  i  te  (6,  22)  f  scntr 

format  < /,  t22,  'c  entral  fs  =',f6.  3) 

establish  the  louer  bound  of  the  capacity-demand  functional 

phimin  =  bound  <  1 ,  1  ) -M-p  i  180 

cmin  =  bound (2, 1 ) 

denmax  =  bound  (3,2) 

ur 1 te  (6, 23) 

format  ( 5 ( / ) , t20,  '  cal  c ulat i on  of  min  fs') 

phil  =  phimin/pilSO 

write  (6,25)  ph i  1 ,  cmin,  denma x 

format  (/,t22,  '  phi  =  ',f6.  2,  /, 

1  t22,  '  cohesion  =  ',  f8.  1,  /, 

2  t22,  '  density  =  ',  f6.  1  ) 

call  fsafty  (slope, hight, phimin, cmin, denmax, fsmin) 

wr 1 te  (6, 26)  f smin 

format  (t22,  'min  fs  =',f6.  3) 

if  (fsr-nin.ge.  1.0)  then 

print*,  '  probability  of  failure  =  0' 

stop 
sndif 
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40 


45 


establish    the       'jpper    bound    of    the    cap<=ic  i  ty-demand 

functional 

phimax  =  bound < 1, 2)«p i 180 

cma  x  =  b  ound (2*  2) 

denmin  =  bound(3»  1) 

lijrite  (6,27) 

format  (5(/).t20,  'calculation  of  rax  f s  '  > 

phi2  =  phimax/pilSO 

ijjT-ite     (6<25)    ph  i2;  cmax,  denmin 

call    fsafty     (slope,  hight,  phimax,  CiTax,  den rTiin.fs max) 

ujTite     (6..  28)     f smax 

format     (t22,   'max     fs    =',  f6.  3) 

if     (f smax.  It.  1.  0)    then 

print*,   '    probability    of    failure    =    1.0' 

stop 
end  1  f 

use  the  rosenbleuth  approximation  to  approxmimate 
the  mean  value  and  the  variance  of  the  capacity- 
demand  functional. 

the  following  expressions  for  the  point  estimates 
of  the  capacity-demand  functional  assume  that  the 
coefficient  of  skeuiness  of  the  phi,C/and  gamma 
parameters  are  all  =  0  and  that  phi,c  and  gamma 
are  unc orre 1 lated 


1 1 


(  10 
WT 

fo 
i  c 
end  if 
s  u  m  f  s 
smf sqr 
do  50 


ut.  eq.  1  )  then 

ite  (6,40) 

rmat  (lhl,t30,  'rosenbleuth  point  estimates') 

ount  =  0 

=  0.  0 

=  0.  0 
i  =  1,  2 


do  50 


do  50 


if  (i 
if  (i 
J  =  1 
i^  (J 
!<=  (J 
k  =  1 
if  (k 


if 
if 


(  k 


eq.  1  ) 
eq.  2) 


eq.  1  ) 
eq.  2) 
2 

eq.  1  ) 
eq.  2) 
(  i  out.  eq. 
icount  = 


phi  = 
phi  = 

cohes  = 
cohes  = 

gamma  = 

gamma  = 

1)  then 

icount 


( mean ( 1  ) 
(mean ( 1 ) 

=  mean (2) 
=  mean ( 2) 

=  mean(3) 
=  mean(3) 

+  1 


stddevd  )  )*pil80 
•  5tddev( 1 ) )*p 1 180 

■  stddev(2) 

■  5tddev(2) 

■  5tddev(3) 

■  stddev(3) 


write  (6,45)  icount 

format  ( 5( /), t20,  'point  est imate ' ,  i3, / ) 
phi3  =  phi/piiao 
write  (6,25)  ph  i3,  c  ohes,  gamma 
end  if 

calculate  the  point  estimator  of  the  factor  of  safety 
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46 


60 


10 
20 

c 
c 


1 


call  fsafty  < si ope< h i gh t- ph i; c oh es . gamma/ fs ) 
uri te  <6>  46)  f s 

format  < t20.  'factor  of  saf ety=  '  <  f 6.  3 ) 
sumfs  =  5umf5  +  fs 
smfsqr  =  smfsqr  +  fs**2 
continue 


calculate  the  mean  value  of  the  point  estimators  of  the  fs 

f smean  =  sumf s/S.  0 

calculate  the  variance  of  the  point  esti.Taters  of  the  fs 

varf s  =  smf sqr/B. 0  -  f 5mean**2 
Sigma  =  sqrt<varfs) 

calculate  the  probability  that  the  fs  u/ill  be  <    f sx  =  l.  0,  i.  e.  ^ 

the  probability  that  the  slope  can  not  maintain  limit  equilibriun 

call  beta  ( f smean.  sigmai  f smin^  f smaxi  f s x , p f a i  1  ) 

write  (6,60)  f  5cntr»  f  smax.  f  smin,  v,  p  f  ai  1 

format  (  Ih  1,  20(  /  ) ,  t20, 

'properties  of  the  capacity-demand  f unc t i onal ',  /i 
:         t30,  'mean  value  ='.  f6.  3,  /, 
3         t30,  'max   value  =',  f6.  3,   /, 
I-         t30.  'min   value  =',f6.  3>/i 
}  t30,  'coef.  of  vav.     =',f6.  3,  /. 

>         t30,  'prob.  of  failure='j  f7.  5) 

stop 

end 


subroutine  fsafty  < slope, h i gh t/ ph i / cohes, gamma, fs ) 

common/hate/  numx, numy, toe, dangle, nsmax 

common/lost/  coti,5inphi 

real  nsmax 

integer  toe, cycles 

data  toll,  tol2/  cycles,  iout/1.  Oe-B,  0.  01,  15,  1/ 


pi  =  acos  (-1.  0) 

if  (lout.  eq.  1)  iijrite(6,  10) 

format  (///,tlO,  'fs  phi',t30,  'fs  cohes',t50,  'stability  number',// 

format  (tlO,  f5.  3,  t30,  f5.  3,  t53,  f7.  4) 


ii=  (ab5(slope-pi/2.  0).lt.  toll)  then 

coti  =  0.  O 
else 

coti  =  1. 0/tan ( slop e ) 
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end  i  f 

smphi  =  sin(phi) 

cgamh  =  c oh es/ < gamma*h i g h t ) 

f  sphi  =  1.0 

f  3C0h3  =  0.  0 

fclold  =  0.  0 

fplold  =  0.  0 


dolOOi=l»cycles 
c 

c      estimate  phi  required  for  limit  equilibrium 
c 

phireq  =  atan < tan ( ph i  ) /fsph i  ) 
sinphi  =  sin<phireq) 
c 

fc2old  =  fclold 
fp2old  =  fplold 
fclold  =  fscohs 
fplold  =  fsphi 
c 

call  phicir 
c 

c      back  calculate  the  fs  on  the  cohesion 
c 

f=coh3  =  cgamh/nsmax 
fratio  =  ab 5 ( f sc oh 5-f sph i ) /f sc oh  5 
c 

c      compare  the  fs  assumed  on  phi  to  the  fs  calculated  on  the  cohesio 
c 

if  (fratio.  It.  tol2)  then 

if  (iout.  eq.  1)  write  (6.  20)  fsphi,  fscohs/ns  ma  x 
fs  =  (fscohs  +  f5phi)/2.0 
return 
elseif  ( fratio.  ge.  tol2)  then 

if  (iout.  eq.  1)  write  (6,  20)  fsphi,  fscohs,  nsmax 
if  (phi.  It.  toll  )  then 

fsphi  =  fscohs 
else 

diffl  =  ab5(fp2old-f sphi ) 

diff2  =  abs(f c2old-f scohs) 

if  (dif f  1.  It.  tol2.  and.  diff2.  It.  tol2)  then 

f sph  i  =  fscohs 
else 

fsphi  =  (fsphi  +  fscohs>/2.  0 
end  i  f 
end  i  f 
end  i  f 
if  '^ieq.  cycles)  then 

print-*,  '  convergence  not  obtained' 
stop 
end  1  f 
1 00    continue 
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end 


subroutine  phicir 

common/hate/  numx<  numy»  toe<  dangle/ nsmax 

common/groovy/  nl imi t< dep th f 

c ommon/ los t/  coti.sinphi 

common/love/  en<  x, y, cscxi  cscy*  cotx<  coty<  stbnum 

real  nsma X/ nl imi t 

integer  toe 


y   is  1/2  the  angle  swept  out  by  the  circle  in  question 

y  =  0.  0 

nsma  x  =  0.  0 

do    25    itery    =    l<nuffly 

y  =  y  +  dangle 

coty  =  1.  0/tan<  y ) 

cscy  =  1, 0/sin( y ) 

X  =  0.  0 

do  25  iterx  =  1/ numx 

x=x+dangle 

iF  (toe.eq.  l.and.x.lt.  y)  go  to  25 

cotx  =  1.  0/tan<  x  ) 

c  sc  X  =  1 .  0/sin ( X  ) 

calculate  the  extent  of  the  limit  eqquilibrium  surface  beyond  the 
toe 

en  =  0.  5*(cotx  -  coty  -  coti  +  sinph i»c sc x*c sc y ) 

if  (toe.  ne.  3.  and.  en.  gt.  0.  0)  go  to  25 

if  ( toe.  eq.  3.  and.  en.  gt.  nl  imit )  go  to  25 

calculate  the  depth  factor 

d  =  0.  5*< CSC x*c scy  -  cotx»coty  +  1.0) 

if  (toe.  ne.  1.  and.  d.  gt.  depthf.  and.  y.  gt.  x)  go  to  25 

calculate  the  value  of  the  stability  number  ffor  the  angles   x  an 

callstabnm 

choose  the  maximum  value  of  the  stability  number 

nsmax  =  ama x 1 < nsmaxi  stbnum ) 

continue 
return 

end 
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subroutine  stobnm 

CO mm on/love/  en*  Xi  y*  cscx*  cscyi  cotxi  coty.  stbnum 

common/lost/  cpti,sinphi 


use  the  friction  circle  method  to  determine  the  stability  numbeT 
of  a  specified  surface  on  the  slope  question 
equations  are  from    taylor  (1937) 

sec  X  =  1.  O/cos  <  X  ) 

toe  circle 

1  f  (  en.  le.  0.  0)  then 

calculate    h/2d  eq.  9# 

paraml  =  ( 0.  5*c sc x*«2»( y *c scy **2  -  coty)  +  cotx 

1  -  coti  )/(l.  0/3.  0*(1.  0  -2.  0'>coti*«-2)  +  coti 

2  «(cotx  -  coty)  +  cotx»coty) 


calculate 


cot (u ) 


eq.   10# 


pora(n2  =  paraml  *y*sec  x»c  sc  x*c  sc  y«-*2  -  tan(x) 
■J  =   a  tan  (  1 .  0/param2) 


calculate 


s  in ( u-v) 


eq.  11# 


param3  =  paraml*s in < u ) *c sc x*c sc y *s inph i 
uv  =  asin(param3) 
V  =  u  -  uv 

calculate  the  stability  number     eq.  12# 

stbnum  =  ( 0.  5*c sc x**2*( y*c 5cy»#2  -  coty)  +  cotx  - 
1  coti )/(2.  0»cotx/tan<v)  +  2.0) 

slip  surface  beneath  the  toe  of  the  slope 

elseif  (en.  gt.0.0)  then 

eq.   14# 

paraml  =  (  0.  5«-c  sc  x**2*(  y*csc  y**2  -  coty)  +  cotx 

1  coti  -  2.  0»en) /(  1.  0/3.  0-*(l.  0  -  2.  0*coti 

2  **2)  +  coti*(cotx  -  coty)  +  cotx*coty  + 

3  2.  0«-en*«-2  -  2.  0*en«-s  inph  i  *c  sc  x<-c  sc  y  ) 

eq.   15# 


piram2  =  paraml  «-y«-5ec  x*c  sc  x*-c  sc  y*-^2  -  taniX) 
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u  =  atan ( 1 . 0/param2) 

eq   16# 

param3  =  paraml*s  in  (  u  )*c  sc  x-»-c  scy«-s  inph  i 
uv  =  as  in ( param3) 
V  =  u  -  uv 

stability  number     eq.  17# 

stbnum  =  ( 0.  5*c sc x**2*< y*c sc y **2  -  coty)  +  cotx 
1  coti  -  2.  0*en  )  /  (2.  0*cotx/tan  (  v)  +  2.0) 

end  i  f 
return 
end 


subroutine  beta  ( xbar< sx^  xmin< xmax/ x> pf ) 
real  const ( 5) 

subroutine  to  find  the  probability  of  given  points 

by  using  a  beta  distribution  fitted  kno'^jing  the 

mean/  standard  deviation  and  range  of  all  data. 

formulas  from   'mechanics  of  particulate  media'  by  m.  e.  harr 

const! 1)  =  xmin 
const(2)=xmax 

calculate   x-tilda    pg.  496   eq.  c-30 

xt=(xbar-const(l) )/(const(2)-con5t(l) ) 

calculate   var iance-ti 1  da  eq.  c-30 

vt=(sx/(const<2)-const(l)) ) **2 

calculate    alpha         eq.  c-31a 

const(3)  =  xt»»2»(l.  0-xt)/vt-(l.  0+xt) 

calculate   beta  eq.  c-31b 

const(4)  =  (con5t<3)+l.  0) / x t-( c onst ( 3) +2.  0) 

const  (5)  =  1.0 

cal  =  qudrtr(con5t(l). const(2),  const) 

c  onst ( 5)=1.  0/cal 

of  =  qudr tr ( const ( 1 )<  X , c onst ) 

return 

snc 
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function  qudrtr ( a<  b/  coef f ) 
real  x g ( 8 ) , wg < 8 ) ,  c oef f ( 5 ) 
data  X  g / ug 

1  /O  0950125098,  0.  2816035507,  0.  4580167/76,  0.  6178762:444, 

2  0  7554044003,  0.  8656312023,  0.  94';5750230,  0.  989400^'349, 

3  0.  1894506104,0.  1826034150,0.  1691565193,0.  1495959888, 

4  0.  1246289712,  0.  0951585116,  0.  062='535239,  0.  0271524594/ 
c 

c 

c      integrates  the  function  betaf  betu-'een  the  limits  a  and  b 

c      by  sixteen  point  symmetric  gaussian  quadrature. 

c      coeff  is  a  real  array  of  coefficients 

c      for  use  in  the  function  betaf,  assumed  independent  of  the 

c      variable  over  which  the  integration  is  taking  place. 

r 
C 

sum=0.  0 

amb  =  (b-a)*0.  5 

apb  =  (b+a)»0.  5 

do  1  i  =  l,8 

xp=apb+amb*xg ( i ) 
xm=apb  — xg  <  i )«amb 

sum=sum+»ug  (i)«<betaf(xp,  coeff)+betaf(xm,  coeff)  ) 
1      continue 

sum=sum*amb 

qudr tr=5um 

return 

end 
c 
c 
c 

function  betaf ( x , coeff ) 

real  coef f ( 5),  x 
c 

c      computes  beta  function 
c  coeff(l)=a 

c  coeff(2)=b 

c  coef f (3>=alpha 

c  coef f (4)=beta 

c  c oef f ( 5)=normal i z ing  constant  at  the  point  x, 

c  where   a  '-I  X  <  b 

c 
c 

c      eq.  c-27a 
c 

betaf  =  coef f <5)»( x-coef f ( 1 ) )**coef f (3)*<coef f (2)-x  )**coef f (4) 

return 

end 
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APPENDIX  F 
STRESSES  CAUSED  BY  CONSTRUCTION  OF  AN  EMBANKMENT 

the  follouing  expressions  may  be  used  to  calculate  the 
stress  changes  caused  by  a  semi-embankment  loading  in  a 
semi-infinite.  u>eightless<  linear  elastic  half-space  (Jut — 
genson,  1940): 


i.c      =-Cj8  +  -a--^  <x-b>>  <F.1> 

z    IT       a       ^ 

^2 


^2 


^2 

ujhere 

P  =  VH 

V  =  embankment  soils  density 

H  =  embankment  height 

B,    a,    a^  b^  r  .  r,,  r^  are  illustrated  in  Figure  F.l. 


Equations  F.l  to  F.3  assume: 

1>   The  ground  surface  is  horizontal. 
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2)  The  embankment  is  perfectly  flexible.   This  is 
equivalent  to  saying  that  the  embankment  - 
halfspace  interface  is  frictionless. 

3)  Strains  are    infinitesimal. 

4)  The  embankment  is  infinite  in  its  linear 
direction. 

The  stress  changes  for  an  entire  embankment  may  be 
obtained  by  summing  the  stress  contributions  of  each  half  of 
the  embankment  with  equations  F.  1  to  F.  3.   It  should  be 
noted  that  the  stresses  cannot  be  calculated  at  positions  A> 
B.  or  C  in  Figure  F.  1. 

Given  ^u    ,    &a    ,    and  At   /  the  principal  stress  changes, 
icr.  and  d>tT^   are   calculated.   The  stresses  in  the  linear 
direction  depend  on  the  values  of  elastic  constants.   Since 
the  problem  is  one  of  plane  strain: 

uihere  V   =  Poisson's  ratio. 

The  value  of  Poisson's  ratio  does  not  affect  the  values  of 

Instructions  for  a  computer  program  that  facilitates  compu- 
tation of  equations  F.  1  to  F.  4  are    included  on  the  follouing 
pages. 
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User  Manual  -  Stress  Under  an  Enbankment 

The  following  program  contained  hereafter  was  developed 
on  a  CDC  6000  series  computer  using  77  FORTRAN.   Input  is 
intended  to  be  in  English  units*  but  other  d imensional ly 
homogeneous  units  may  be  employed.   All  input  is  unformat- 
ted. 

1.   Read  in  on  one  record: 

a)  the  side  slope  of  the  right  hand  side  of 
the  embankment  (degrees) 

b)  the  side  slope  of  the  left  hand  side  of 
the  embankment  (degrees) 

c)  the  embankment  height  (ft) 

d)  the  crest  to  crest  width  (ft) 

e)  the  density  of  the  embankment  soil  (pcf) 

2    Read  in  on  one  record: 

a)  the  Poisson's  ratio  of  the  soil 

3.  Read  in  on  one  record: 

a)  the  number  of  ( Xi 2 )  c oord inate  pairs  at 
which  stresses  are  desired  (integer) 

4.  Read  in  on  one  record: 

a)  the  X  coordinate  of  the  point  at  which  stresses 
are  desired  ( f t ) 

b)  the  Z  coordinate  of  the  point  at  which  stresses 
are  desired  ( ft ) 
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Repeat  for  each  coordinate  pair. 

X  is  measured  from  the  center  line  of  the  crest 
of  the  embankment  (see  Figure  F.  2). 

Z  is  measured  dountusrds  from  the  ground  surface 
( see  Figure  F.  2). 


The  program  listing  is  included  in  the  follouing  pages. 


206 


-^N 


< 

o 
o 
o: 

Q. 

CO 
CO 
LJ 

cr 
co 

ir 
o 
u. 

>- 
q: 

I- 

UJ 

o 
o 

H- 
Z 
IxJ 


< 

m 


CVJ 


Ul 

q: 


207 


program  main  < inputs  output/  tape5=input,  tap e6=outp ut ) 

common  x,  z,  hembnk<  uiidth,  densty.  sigmaz,  sigmax,  sigmxi,  sigmal. 
1        5igma2>  sigma3,  thetaimui  slopel.  5lope2,  pi 


real  mu 


this  program  calculates  the  values  of  stress  change  caused 
by  an  embankment  on  a  linear  elastic  half-space. 


read  in 

1)  the  righthand  side  side-slope  of  the  embankment  in 
degrees. 

2)  the  lefthand  side  side-sloe  of  the  embankment  in 
degrees. 

3)  the  embankment  hight 

4)  the  crest  uiidth 

5)  the  embankment's  density 

read  (5,*)  slope  1 ,  si ope2/  hembn k . uid th.  d ensty 

write  (6<5)  slope  1*  5lope2,  hembnk ;  ji d th ,  d ensty 

format  (  '1  ',  5</),  t50.  'right  side-slope  =  '<f8.  2,  '  degrees.  ' ,  / , 

1  t50.  'left  side-slope  =  ',f8.  2,  '  degrees.  ',/, 

2  t50,  'embankment  hight  =  ',flO.  2-  '  feet.  ',/, 

3  t50,  'crest  width  =  ',  f6.  2>  '  feet.  ',  / , 

4  t50,  'embankment  density  =  ',f7.  2.  '  lb./ft**3') 
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pi  =  acos(-l.O) 

read  in  the  value  of  poissons'  ratio  for  the  elastic 
hal f-space. 

read  (  5i  *  >  mu 

tur  1  te  (  6<  3  )  mu 

format  (//,  t50,  'mu  =  ',  f4.  2) 

read  in  the  number  of  coordinate  pairs  that  the  stress 
change  is  desired  at. 

read  (  5.  »  )  npairs 

ur  i te  (6.  10) 

-format  (  '  1  ',  tS.  '  x  S  t22/  '  z' ,   t38.  'sigma  z     ',  t52,  '  sigma  x  '-  t66, 

1  'Sigma  xz  '»  tSl/  'sigmal  '  >  t95,  'sisma2'/  tl09,  'SigmaS'. 

2  tl24>  'theta  ', //) 


calculation  loop    input  an  (x,z)  coordinate  pair 
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for  each  loop 


do  100  i  =1/ npairs 

read  (5,  »)  x,  z 
call  stress 

urite  <6/20)  x,  Z/  s  i  gmaz  <  si  giTa  x  /  si  g.Tix  z.  s  i  gma  1  -  s  i  gma2-  s  i  gmaG 
1  theta 

20  format  (2  (  f  9.  1,  6x  ) ,  6  (  5x,  f  9.  1  > ,  5x  ,  f  9  2) 

100    continue 
end 

subroutine  stress 

common  x-  z, hembnk, width,  den  sty.  sigmaz,  sigmax,  sigmxz-  sigmal, 
1        sigma2>  sigma3. theta. mu.  slopel.  5lope2.  pi 


real  mu 


p  =  densty  ♦  hembnk 

c 

c 

c  contribution  to  the  stresses  from  the  right  half  of  the 

c  embankment  loading 

c 

c 

sloper    =    slopel    *    pi/180.0 

ar    =    hembnk/tan ( sloper ) 

br    =    ar    +    ujidth/2.  O 

X 1    =    br    -    X 

r2    =    sqrt  (  (ijidth/2.  0   -^   ar    -    x  1  )  **2    +    z«*2) 

rl    =    sqrt((xl    -    ar)*»2    +    z»*2) 

rO    =    sqrt(xl**2    +    z»«2) 

betar  =  acos<<rl**2  +  r2**2  -  ( wi d th /2. 0) **2 ) / (2. 0*r l*r2 ) ) 

gammar  =  acos((rl**2  +  r0»*2  -  ar --*2  )/(  2.  0*r  l*rO  )  ) 

ciflr  =  betar  +  xl/ar  *  gammar  -  z* ( x 1-br ) /r2**2 

rightx  =  betar  +  x l/ar*gammar  +  z»( x 1-br ) /r2**2  +  2. 0*z/ar»al og 
1  (rl/rO) 

rihtxz  =  z/ar*gammar  -  z**2/r2*»2 
c 
c 

c  contribution  to  the  stresses  from  the  left  side. 

c 
c  - 

slopel  =  slope2*pi/180.  0 

al  =  hembnk/tan ( slopel ) 

bl  =  al  +  width/2.  0 

x2  =  bl  +  X 

r2  =  sqrt(  (ujidth/2.  0  +  al  -  x2)»*2  +  z**2) 

rl  =  sqrt((x2  -  al)»*2  +  z**2) 

rO  =  sqrt(x2-R-*2  +  z**2) 

betal  =  acos ( <rl**2  +r2**2  -  ( wi d th /2. 0 ) **2 ) / (2. 0*r l*r2 > ) 

gamma  1  =  ac  os  <<  r  1**2  +  r0*#2  -  al**2;' /(  2.  0-«-r  l*rO  .''  ) 

cifll  =  betal  +  x2*gammal/al  -  z*  (  x2-b  1  >  .'■r2*-*2 

leftx  =  betal  +  x2/al*gammal  +  z/r2*»2*<v2  -  bl)  +  2.  0*z /al*al oc 
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1  (rl/rO) 

leftxz  =  z/al-K-gammal  -  2**2/r2**2 
sigmaz    =  p/pi»(ciflr  +  cifll) 


sigmax  =  p/p i »(rightx  +  leftx) 
sigmxz  =  p/pi»(rihtxz  +  leftxz) 
descr  =  sqrt  ((sigmax  -  si gmaz )**2/4. 

=  (sigmax  +  5igma2)/2.  0 

=  center  ■*•    descr 

=  center  -  descr 

=  mu*(5igmal  +  sigma3) 


s  i  gmx  z**2) 


center 
s  i  gmal 
s 1 gma3 
si  gma~ 


theta  =  (0.  5»atan(-2.  0»sigmxz/( sigmaz  -  s igmax ) ) )*180.  0/p i 


return 
end 
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APPENDIX  G 
MAGNITUDE  OF  CONSOLIDATION  SETTLEMENT 

The  program  contained  hereafter  coi<tputes  the  magnitude 
of  consolidation  settlement  of  compressible  soil  layers 
caused  b\i    construction  of  an  embankment.   Up  to  ten  layers 
are    permitted.  Each  layer  may  be  automatically  divided  into 
any  number  of  strata.  The  preconsol idation  pressure  profile 
is  input  as  shoun  in  Figure  4. 2b.   The  void  ratio  is 
automatically  corrected  for  depth  uith  equations  4.6a  to 
4.  6h.   Settlement  is  computed  luith  equations  4.2  to  4.5. 

The  program  was  developed  on  a  CDC  6000  series  computer 
using  77  FORTRAN. 
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User  Manual  -  Magnitude  of  Consolidation  Settlement  Program 

Input  is  intended  to  be  in  English  unitsi  but  other 
d  iniensional  ly  homogeneous  units  may  be  employed. 

All  input  is  unformatted. 


Specify  the  Number  and  Type  of  Layers  in  the  Soil  Profile: 

1.  Read  in  on  one  record: 

a)  the  total  number  of  soil  layers  (integer) 

b)  the  number  of  compressible  soil  layers  (integer) 

c)  the  number  of  layers  with  a  portion  above 
the  ground  uater  table  (integer) 

d)  the  number  of  layers  uiith  a  portion  below 
the  ground  tuater  table  (integer) 

Specify  the  Compressibility  Models  for  Each  Layer: 

2.  For  each  soil  layer,  read  in  on  one  record: 

a)  the  layer  number  (integer) 

b)  the  thickness  of  the  layer  (feet) 

c)  An  identifier  variable  that  indicates 

if  the  soil  layer  is  compressible  (integer). 

IF  the  soil  is  not  considered  to  be  compressib lei 

input  '0' 

If  the  soil  is  considered  to  be  compressible. 

input  '  1  '. 

d)  If  the  soil  is  considered  to  be  compressible/  read 


212 


in  the  follouiing  variables  on  a  separate  record 

prior  to  inputting  (a)-(c)  for  the  subsequent  soil 

layer. 

i)  An  option  variable  that  is  used  to  control  the 

manner  by  uihich  e   and  o"      are  assumed 

to  vary  with  depth  in  the  soil  layer  (integer) 

If  one  value  of  e   and  o"      ar-e   used  to  repre- 
o        p 

sent  the  entire  soil  layer,  input  '1'. 

If  e   and  o'  &Te   specified  at  three  depths 
o        p       "^ 

within  the  soil  layer<  input  '2'. 

ii)  An  option  variable  that  specifies  if  the  soil 

is  to  be  treated  as  normally  consolidated  or 

underconsol idated  at  depths  beneath  which  the 

calculated  value  of  o"^     is  greater  than 

the  value  of  ct'   <inteQer> 
P 

If  the  soil  layer  is  to  be  treated  like  normally 

consolidated  soil  at  depths  where  a'   is 

greater  than  o"    ,    input  "1^. 
P 

If  the  soil  layer  is  to  be  treated  as  an  under- 
consolidated  soil  at  depths  where  o'   is  greater 

than  the  value  of  a'  ,    input  '2' . 

P 


Specify  the  Ground  Water  Table: 

3.   Read  in  on  one  record: 

a)  the  depth  of  the  groundwater  table  beneath 
the  ground  surface  (feet). 
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The  depth  is  measured  douinuards  from  the 
ground  surface  (see  Figure  F. 2). 

Spscify  the  Saturated  Density  of  the  Soil  Lagers: 

4.  For  each  lager  with  a  saturated  zone/  read  in  on  one  record: 

a)  the  lager  number  (integer) 

b)  the  saturated  density  of  the  soil  layer  (pcf) 

Specify  the  Unsaturated  Density  of  the  Soil  Layers: 

5.  For  each  layer  with  an  unsaturated  zor.ei  read  in  on  record: 

a)  the  layer  number  (integer) 

b)  thedensity  of  the  soil  layer  (pcf) 

Specify  the  Compressibility  Parameters  of  Each  Layer: 

6.  For  each  compressible  layer<  read  in  on  one  record: 
a)  the  layer  number  (integer) 


If  one  value  of  e   and  o""   represent  the  entire  layer 
op 

read  in  on  a  separate  record: 


b)  the  initial  void  ratio 

c)  the  preconsol idation  pressure  (psf) 

d)  the  compression  index 

e)  the  recompression  index 

f)  the  depth  beneath  the  ground  surface  from 
which  the  tested  sample  was  taken  (ft) 
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If  the  soil  may  be  considered  to  be  undercotisol  idated  /  input 
on  a  separate  record: 

b>  the  compression  index 
c)  the  recompression  index 

Also/  if  the  soil  is  underconsol i dated,  input  the  following  irems 
on  the  subsequent  record: 


d>  e   at  depth  1# 
o 

e>  n'       at  depth  1#  <psf> 
f)  Depth  1#  (ft) 


g!-"  e   at  depth  2# 

h>  o"       at  depth  2#  Cpsf? 
P 

i)  Depth  2#  (ft) 


j  >  e   at  dept  h  3# 

k>  Q"       at  depth  3#  <psf> 

1)  Depth  3#  (ft) 


Specify  the  Embankment  Load: 


Read  in  on  one  record: 

a)  the  height  of  the  embankment  (feet) 

b)  the  sideslope  angle  of  the  embankment 
(degrees  ) 

c)  the  uidth  of  the  embankment  crest  (feet) 

d)  the  density  of  the  embankment  soil 
(pcf  ) 


215 


Specify  the  Lateral  Limits  of  the  Settlerrent 
Calculations : 

8.  Read  in  on  one  record: 

a)  the  leftmost  bound  for  which  settlements 
will  be  calculated  <feet) 

b)  the  rightmost  bound  for  which  settlements 

will  be  calculated  (feet) 

c)  the  number  of  evenly  spaced  points  along  the 

embankments  cross  section  at  which  settlement 
will  be  calculated  (integer) 

Coordinates  for  these  bounds  are  measured  (+)  and  (-) 
to  the  right  and  left  of  the  embankrnent  centerline 
respectively. 

If  the  settlement  of  only  one  profile  is  desired/  the 
program  calculates  settlement  at  the  left  boundary. 

Specify  the  Number  of  Strata  Compressible  Layers  are 
Divided  into: 

9.  For  each  compressible  layeri  read  in  on  one  record: 

a)  the  layer  number  (integer) 

b)  the  number  of  strata  into  which  the  layer  will  be 
divided  into  for  purposes  of  analysis  (integer) 


The  program  listing  is  contained  in  the  following  pages. 
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parameter  (nprfil  =  11-  nlyrs  =  10,  nlyrsl  =  11) 

integer  compr5<  nlyrs  ),  5trata(nlyrs)i  ocr(nlyrs).  iprfiUnlyrs), 
1         i  out 

real  thick(nlyrsl)<  densatCnlyrs).  denmst(nlyrs), 

1  cc  (nlyrs),  totaKnprfil),  settle  (nprfil,  nlyrs), 

2  re (nlyrs),  xll  (nprfil),  eol (nlyrs),  eo2(nlyrs), 

3  eo3(nlyr5), pcl(nlyrs), pc2(nlyr5), pc3( nlyrs),  zsampl(nlyrs) 

4  2samp2(nlyrs )/  zsamp3(nl yrs ) 

data  gammaiu/62.  4/,  tol/.  001/ 

program  assumes  e-logp  behaviour  only  in  saturated  zone. 


units  of  length  -  (feet) 
units  of  force   -  (pounds) 
units  of  pressure  -  (psf) 


Lur  i  te  (6,  1  > 

wr  i  te  (  6,  2  ) 


pi  =  acos(-l.O) 

read  in: 

1)  total  #  of  soil  layers 

2)  #  of  compressible  soil  layers 

3)  #  of  layers  with  portions  above  the  ground  niater  table. 

4)  4*  of  layers  uith  portions  beloin  the  ground  water  table. 


read  (5,«)  nlay er, nc layr , nabove, nbe 1 ou 


do  4  i=l, nlayer 

ocr ( i  )  =  0 

iprf il  (  i  )  =  0 
continue 

read  in  geometry  of  problem 


for  each  layer  read  m: 

1 )  lay  er  # 

2,'  the  thickness  of  the  I'th  layi 
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3)  if  the  layer  is  compressible, 
if) yes  then  input   '1' 
if, no   then  input   '0' 
&)     if  the  layer  is  compressible,  read  in  on  a  seperate 
card  : 
iprf il 

if  the  variation  of  pc  and  eo  is  defined  by 
two  line  segments,   set  iprfil  =  2 
if  only  one  value  of  pc  and  eo  are    specified 
for  the  layer,  set  iprfil  =  1 
ocr 

if  pc  is  to  be  auto.ratically  corrected  to 
the  overburden  pressure  ujhen  pc'Iover  b  ur  den 
pressure,   set  ocr  =  1 

if  pc  is  to  be  considered  urd ere onso  1  i dated 
luhen  pc  <  overburden  pressure, 
set  ocr  =  2 

ocr  =1       if  the  soil  is  normally  consolidated 
at  depths  beneath  uihich 

the  overburden  pressure  is  greater  than  the 

value  of  the  prec onsol i dat i on  pressure  thai 

luas  input 
ocr  =  2      if  the  prec onso  1  i dt ion  pressure  is  constan- 

with  depth. 


do  5  i=l,nlayer 


10 


read  (5,«)  i  1 ,  th  i  c  k  (  i  1  ) ,  c  omprs  (  i  1  ) 


if  (  comprs  (  i  1  ) .  eq.  1  )  read  <5,*)  i  pr  f  1 1  (  i  1  )  ,  ocr  (  i  1  ) 


continue 

do  10  i=l,nlayer 

write  (6/6)  i,thick(i) 

if  (comprs(i).  eq.  1)  then 

if  (ocr(i).eq.  l.or.  iprfil(i).eq.  2)  then 

wr  ite  (6,7) 
elseif  ( ocr ( i  ) .  eq.  2)  then 

write  (6,8) 
endif 
elseif  <comprE(i).eq.  0)  then 

write  (6,9) 
end  i  f 
continue 
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read  in  elevation  of  ground  water  surface 


read  ( 5, *)  zwater 
ur  i te  (6j  1 1 )  zwater 


24 


ijjr  ite  (6.  20) 


initialize  moist  and  saturated  densities  of  each  layer. 


do  24  1=1/ nlayer 

denmst  (  i  )  =  0.  0 
densat  <  i  )  =  0.  0 

continue 


c 
c 
c 
c 
c 
c 
c 
c 

c 
c 

c 
c 
25 

c 
c 
c 
c 
c 
c 
c 

c 
c 


for  each  layer  with  a    saturated  zone-  read 

1 )  layer  # 

2)  saturated  density 


do  25  i  =  l  /  nbel  ou) 


read  (5j«)  i  1>  densat  (  i  1  ) 


continue 


for  each  layer  with  an  unsaturated  lonei  read  in: 

1)  layer  « 

2)  moist  density 


do  30  i=l ; nabove 


read  <5<*)  i 1 . denmst ( i 1 ) 
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c 
30 

c 
c 
c 
c 
c 
c 

c 
c 
c 
c 
c 

c 
c 
c 

c 
c 
c 
c 
c 
c 
c 
c 
c 

c 
c 
c 

c 
c 


c 
c 

40 

c 
c 
c 
c 
c 


continue 


read  in  compressibility  properties 


do  40  i  =  1, nclayr 


input  the  number  of  the  compressible  layer 


read  (  5.  *  )  i  1 


if 


(  iprf  11  (  i  1 ).  eq.  1 )  then 

for  each  compressible  layerj 


read  in: 


1)  initial  void  ratio 

2)  preconsol idation  pressure  (psf) 

3)  compression  index  (log  10) 

4)  recompression  index  (log  10) 

5)  sample  depth  (feet) 


read  (5,*)  eol  (  i  1)  i  pc  1  (  i  1  ) ,  c  c  (  i  1  )  ,  r  c  <  i  1  )  <  z  samp  1  (  i  1  ) 


elseif  (  iprf  i  1  (  i  1  ).  eq.  2)  then 


read  (5. »)  cc  (  i 1  ) , re ( i 1 ) 

read  (5,*)  eol  (  i  1) ,  pc  1  (  i  1  ),  z  samp  1  (  i  1  ),  eo2  (  1 1  ),  p  c2  (  i  1  ) . 
zsamp2( il)/ eo3(il)<  pc3(il);  zsamp3( i 1 ) 


end  i  f 
continue 


output  material  parameters 


write  (6, 50) 

dc  60  i  =  l,nlayer 

if  (abs(denmst  (  i  )  ).  It.  tol  )  then 

urite  (6(54)  i/den5at(i) 
elseif  (ab5( densat ( i ; ).  It.  tol )  then 
u;rite  (6.56)  i;denmst(i> 


60 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
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elseif<densat<i).gt.O.  0.  and.de nni5t(i:'  gt.O.O)  -then 
write  <6;57)  i  /  d  ensat  (  i  )  <  d  en/nBt  <  i  ) 

mpr 5  (  i  ) .  eq.  1  )  then 

iprf  i  1  <  i  ) .  eq.  1  )  then 

ite  (6.  58)  eol(i)/pcl(i)<cc(i;'<rc(i),zsampl(i) 

eq.  2)  then 

eol(i)<pcl(i)<cc(i),rc(i)/Zsampl(i), 

eo2(i)<pc2(i))  zsamp2<i), eo3(i). pc3(i) 

z  samp3( i ) 


end 

if 

if 

(CO: 

i 

f  < 

elseif 

ur 

end 

if 

e  Ise  i 

f  ( 

end 

If 

c  on 

tinue 

hJT  i 

te  (6, 

67) 

(iprf il(i ) 
ite  (6,61) 


comprs(i).  eq.  0)  then 
write  (6,  59) 


read  in  dimensions  of  embanknent  load: 

1)  embankment  hight 

2)  sideslope  of  embankment  (degrees) 

3)  crest  width 

4)  density  of  embankment  material  duririo  period  of 
consolidati  on. 


read  (5,*)  hembnk»  si  ope  1 , wi d th>  densty 
lurite  (6/68)  hembnk  ,  s  1  op  e  1  /  wi  d  th  >  d  enstg 


read  in : 

1)  the  lateral  limits  between  uhich  settlement  is  calculated 

2)  the  #  of  points  between  the  bounds  for  which 
consolidation  is  to  be  calculated. 


read  (5<-»)  b  ound  Ij  boundr ,  nntrvl 


nintrv  =  nntrvl  —  1 

write  (6j72>  nntrvl,  bound  1 /  boundr 

if  (nntrvl.  gt.  l)then 

d x  =  (boundr  -  b ound 1 ) /f 1  oat  (nintrv) 
else  if  (nntrvl.  eq.  1)  then 
dx  =••  0.  O 

5  n  d  i  f 

wr 1 te  ( 6. 74 ) 
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c 
c 
c 

c 
c 
c 
c 
c 

c 

c 

c 
c 

so 

c 
c 
c 
c 
c 
c 
c 

c 
c 
c 


for  each  compressible  layer)  read  in: 
1  )  layer  # 

2)  the  #  of  strata  the  layer  :j;ill  be  divided  into  for 
purposes  of  analysis 


do  SO  i=l/ nc layr 


read  <5,*)  i  1.  strata  (  i  1  ) 


write  <6<75)  i  li  strata  (  i  1  ) 


continue 


if  strains  are  to  be  output  for  each  strata 

set  iout  =  1 

otherwise   set  iout  =  0 

read  (  5/  ■«■ )  i  out 


X 1  =  b  ound 1  -  d  x 


if  (iout.  eq.  1)  lurite  (6<85) 

this  loop  sums  settlement  from  each  layer   beneath 
each  point  for  uhich  settlement  is  desired. 


do  200  np=  l,nntrvl 

2  top    =0.  0 

pr  top    =    0.0 

xl=xl+dx 

if    (iout.eq.  1)    write    <6i90)     xl 

X    =    abs     ( xl ) 

total  (np  )    =    0.  0 


this  loop  calculattes  settleT.ent  of  each  layer 

beneath  the  np'th  layer, 
do  190  i=l/nlayer 

if  (iout.eq.l)  write  (6,91)  i>thick(i) 
settle  (np.i)  =  0.0 


zbottm  =  ztop  +  thick(i) 
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c 

c  calculate  the  pressure  at  the  bottoiT  of  the 

c  i'thlayer. 

c 

if  (z'jjater  gt.  ztop.  and.  ziuater.  It.  zbottm)  then 

prbotm  =  prtop  +  ( zwater-z top )*denm5 t t i /  +  ( z b ottm-zwater ) 
1  •«•(  densat  <  i  ) -gammau) 

else  if  (  zujater.  le.  ztop)  then 

prbotm  =  prtop  +  (  z  bottm— z  top  )  *  (  d  ensat  (  i  ) -gammaui  > 
else  if  (zuater.  ge.  zbottm)  then 

prbotm  =  prtop  +  (  zbottm-z  top  ) -*c!  enms  t  (  i  ) 
endif  I 

I 

if  <comprs(i).  eq.  1)  then 

dz  =  ( zbottm-z top ) /fl oat ( strata ( i  )  ) 
c 

if  (  iprf  i  1  <  i  ) .  eq.  1  )  then 
c 

c  calculate  the  sample  pressure 

c 

if  (zujater.  le.  ztop)  then 

psampl  =  prtop  +  (densatCi)  -  gamma^ij  )*(  z  samp  1  (  i  ) -ztop  > 

else  if  (zwater.  gt.  ztop.and.zwater.  It.  Z5arT.pl(i))  then 

psampl  =  prtop  +  denmst  (  i  )  *<  zii»ater-z  top  )  +  (  denmst  (  i  ) -gammau;)  * 
1  (  zsamp  1  (  i  )-zu)ater  ) 

end  i  f 

erid  1  f 
c 

nstrta  =  strata(i) 
c 
c 
c 

c  this  loop  divides  each  compressible  layer  into 

c  distinct  strata  for  purposes  cf  calculating 

c  settlements. 

c 

do  100  n=l; nstrta 

z    =    ztop    +    f loat<2»n-l )/2. 0*dz 

if     (iout.  eq.  1)    lurite     <6,92)    n,  z 
c 
c 

c  calculate    the    insitu    pressure    at    the    midpoint  ^ 

c  ofthen'thstrata.  | 

I 
if    ( zwater.  le.  ztop)    then 

prl    =    prtop    +    (  densat  (  i  ) -gamma'j; )  *  (  z-z  top  > 

else    if    (zuater.  gt.  itop.  and.  ziuater.  Itz)     then 

prl    =    prtop    ■*-    (  2Ujater-z  top  ) -f-a  enms  t  (  i  )     +     i.  z -zujater  )  ■»■ 

1  <  densat  <  i  ) -gammauj ) 

endi  f 


if     (.  z .  le.  zuater)    then 
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set    the    settlement    of    the    strata    =    0      above    the    ground 
ii;ater    table. 


r  qee  z  e    =    0.0 


elseif  ( zwater.  It.  z)  then 


calculate  eo  at  the  center  of  each  strata 
beneath  th  mater  table. 


if  (  iprf  i  1  (  i  ) .  eq.  1  )  then 

if  (  ocr ( i ) .  eq.  2)  then 

precon  =  pcl(i) 
elseif  (ocr(i).eq.  l.and.prl.gt.  pcl(i))  then 

precon  =prl 
elseif  (ocr(i).eq.  l.and.prl.le.  pcl(i)>  then 

precon  =  pcl(i) 
end  if 

center  of  strata  above  the  samples  elevation 

if  ^z.  It.  zsampl<i))  then 

if    (psampl.  le.  precon)  then 

case   i 
the  soil  is  over  consolidated  above  the  sample 
estrat  =  eol(i)  +  re ( i ) *al og 10( p sanp 1 / prl ) 


elseif  (psampl.  gt.  precon)  then 
if  (  ocr  <  i  ) .  eq.  2)  then 
if  (prl.  gt.  precon)  then 

case   ii 
the  soil  is  underconsol i dated  and  therefore  uill 
have  a  constant  void  ratio  at  pressures  above 
the  preconsol idat i on  pressure, 
estrat  =  eol  ( i ) 


else  if  (prl.  le.  precon)  then 

case   i 1 i 
tne  soilis  under c onso 1 i dated  at  pressures  above 
the  p'c  and  over c onsol i dated  at  pressures  be- 
1  ouj  the  p'c 

estrat  =  eol(i)  +  re ( i ) *aloe 10( pTecon/pr 1) 
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end  i  f 

elseif  (ocr(i).eq.  1)  then 

case   iv 
the  soil  is  normally  consolidated  at  pressures 
above  the  p'c  and  over-consolidated  at  pressures 
b  eneath  the  p'c 

estrat  =  eol(i)  +  cc ( i )*al og 10( psamp 1 /pr ec on ) 
+   rc(i)*aloglO<precon/prl) 


end  if 
end  if 

center  of  the  strata  belou  the  samples  elevation 

else  if  (z.  ge.  Z5ampl<i))  then 
if  (psampl.  It.  precon)  then 
if  (prl.  le.  precon)  then 

case   V 
the  soil  is  over-consolidated 

estrat  =  eol(i)  -  re ( i ) *a 1 o^ 10 i p r 1 /p samp  1 ; 


else  if  ( prl.  gt.  precon )  then 
if  (  ocr  (  i  ) .  eq.  1  )  then 

case   vi   a ) 
the  soil  is  over-consolidated  beloiu  the  p'c  and  r-ormally 
consolidated  above  the  p'c 

estrat  =  eol(i)  -  r c ( i )*a 1 og lOf precon/psamp 1 )  - 


1  cc ( i )*alog 10( pr 1 /precon ) 

elseif  (  ocr  <  i  ) .  eq.  2)  then 


case  vl  b ) 
the  soil  over-consolidated  belotu  the  p'c  and  over  c  onso  1  i  date 
above  the  p'c 


S5  trat  =  eol  (  i  > 
3nd  1  f 


-  re  (  i  )  «-a  1  og  10  (  prec  on/psamp  1  > 


end  1  f 
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else  if  (  psamp  1.  ge.  precon  )  then 

case   vi  i 
the  soil  is  underc onsol i dated 


estrat  =  eol ( i  ) 


end  if 


end  if 

elseif  (  iprf  il<  i  ).  eq.  2)  then 

if  (  2.  le.  25amp2(  i  )  )  then 

precon  =  pc2(i)  +  <pc2(i>  -  p cl ( i ) ) / < z samp2( i ) 
1  zsamp 1 < i ) )»( 2  -  zsamp2(a)) 

estrat  =  eo2(i)  +  (eo2(i)  -  eo 1 ( i ) ) / ( i samp2 ( i ) 
1  zsampl ( i ) )*( z  -  zsamp2(i>) 

elseif  <  z.  gt.  Z5amp2(  i  )  )  then 

precon  =  pc2(i)  +  (pc3(i)  -  p c2 f i ) ) / ( z 5amp3 ( i ) 
1  Z5amp2(  i  )  )*(  z  -  zsanip2(i)) 

estrat  =  eo2(i)  +  (eo3<i)  -  eo2( i ) ) / ( zsamp3( i ) 
1  Z5amp2< i ) )*< z  -  zsamp2(i)) 

end  if 
if  (  precon.  It.  prl.  and.  ocr  <  i  ).  eq.  1  )  precon  =  prl 
end  if 

call  stress  < x.  z, hembnk,  slope  1, wi d th / densty,  s i gmaz ,  p i ) 
sumsig  =  sigmaz  +  prl 


calculate  the  compression  of  the  strata 


if  (sum 
sqe 

else  if 
if  ( 


else 


sig.  le.  precon)  then 

eze  =  dz/(1.0  +  estrat ) *rc ( i ) *a 1 og 10 ( sums ig /pr 1 ) 

( sumsig.  gt.  precon )  then 
prl.  le.  precon)  then 

sqeeze  =  (dz/<1.0  +  estrat ))*( re ( i ) *a 1 og 10 ( prec on/prl ) 
cc<i)*aloglO(su(7isig/precon)  ) 
if  (prl.  gt.  precon)  then 
if  <  ocr  (  i  ).  eq.  1  )  then 

sqeeze  =  d2/(1.0  +  estrat )*( c c ( i ) *a log 10( sums i g/pr 1 ) ) 
elseif  (  ocr  (  i  ).  eq.  2)  then 
sqeeze  =  dz/(1.0  +  estrat >*( c c ( i > »al og 10 

(sumsig/precon) ) 
endif 


end  i 

endif 


end  if 


if  (  1  out  eq.  i  )  then 

strain  =  sqee z e/d 2*100. 0 
uirite  (6.97)  strain 

end  if 
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c 
c 

c 
100 

c 
c 

c 
c 
c 
c 
c 


c 
c 
c 
c 

190 
200 

c 
c 
c 
c 


300 

c 
c 
c 


400 

c 

c 
c 
c 
c 
c 

1 


sum  the  settlements  of  each  strata  ir  t^le  layer, 
settle  (np.  i)  =  settle  (npi  i)  +  sqeeze 
continue 

end  if 


reset  the  elevation  and  pressure  at  the  top 
of  the  layer  beneath  the  current  layer. 

z top  =  z  b  ottm 
prtop  =  prbotm 

sum  the  settlements  of  each  layer  Der.eath  the 
the  point  at  uihich  the  settlement  is  desired. 

total(np)  =  total(np)  +  settle<np»i) 

c  ont inue 

continue 


output  statements 


write     (6,210) 
xll  (  1  )  =  boundl 

if  (nntrvl.  ge.  2)  then 
do  300  np  =  2* nntrvl 

xll  (np  )  =  xll  (np  -  1 )  +  dx 

continue 
end  i  f 

ujrite     (6-250)     (  x  1 1  (np  ),  np  =  l  <  nntrvl  ) 

lur  i  te    (6;  305) 

do  400  i  =  l<  nlayer 

if  (comprs<i).  eq.  1)  then 

urite  (6»310)  i.  ( sett  1 e (np ,  i )- np  =  l , nntrvl ) 

end  if 
cont  inue 

write  (6/410)  ( total  (np  )*  np  =  1, nntrvl) 
formats 


f  ormat(  '1  ',  25(/)/  t30,  66(  '*'),  /, 
1         t30,  '■»',  t95,  '*',  /,  t30,  '«',  t95,  '*',  /,  t30,  '*',  t95;  '*',  /, 
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2.  t30,  '*'/t55,  'program  elogp'-t??-  '»',/, 

3        t30,  '«', t40,  'this  program  calculates  the  settlement  beneat 

5  t30,  '*'.  t4o!  'an  infinitely  long  e.Tibankment  load  due  to  the 

6  -  t95-  '*'-/.  ,     ,        /  j-o-n 

7  t30,  '«'.  t44,  'compression  oi    saturated  clay  layers.  , tY3, 

9        t30'/»'.t95,  '♦',/,  t30,  '»',t95,  '*',/■  t30,  '»',t95,  '*',/■  t30, 

1        66('*')) 

2      formate  'l'.5</),  t58.  'problem  geometry '.  5(/).  ki=-' 

1  t38,  'layer', 15x,  'thickness  (  f eet )  ' , Sx.  'compressib le  :  - 

2  //) 

6  format  ( t39,  i2»  t60,  f  6.  1) 

7  format  ('+'.  tee.  'yes') 

8  format  ('  +  ',  t32,  'underconsol i dated  '  ) 

9  formate  '+'-  tea,  'nO')  ^  .,,/,, 
20  format(5(/), t52,  'properties  of  the  soil  profile  ./////) 
50  format  ( t2,  '  lay  er ',  tlO.  'sat.  density'- 

1  5x-  'moist  density 'i  7x,  'void  ratio', 

2  5x,  'preconsolidation',  5x,  'compression ',  1  x,  'recompression 

3  ,  5x,  'sample  depth',/, 

4  t49.  'initial '.  t64,  'pressure ',  ta6,  '  inde  x  ' ,  1 100, 

5  'index',/,  tl3,  '(pcf)  ',t31,  '(pcf)  ',  t67,  '  <  psf  )  '-  1 1 16, 

6  '(feet) ', 5(/)) 

5<\  format(t3,  i2,tl3,  f5.  I,t31,  '»»♦»♦') 

56  format(t3,  12,  tl3,  '*»**♦''  t31,  f5.  1) 

57  format(t3,  i2,  tl3,  f5,  1,  t31,  f5.  1)  ^  .=  ,  ^i,^  ^^  in 

53     formate-^',  t50,  f5.  3,  t66,  fS.  1,  t86,  f5.  3,  tlOO,  f5.  3,  Vil;,'  ,!!»** ' 
59     formate  '+',  t50,  '»*»♦*'.  t66,  '«»***>**''  t86,  '»«***'.  tlOO,  »*♦♦*  - 

11     format(3(/>!35x,  'the  phreatic  surface  lies'-f8.  1,'  feet  beneath 

Ithe  ground  surface.  ') 
61     format  (  '+ ' ,  tSO-  f  5.  3.  t66,  f  8.  1 ,  tS6,  f  5.  3,  1 102,  f  5.  3, 

1  tll7.  f6.  1,  /,  t50,  f5.  3,  t66,  f8.  1,  tll7,  f6.  1,  /, 

2  t50,  f5.  3,  t66,  fS.  1-  tll7,  f6.  1> 

67     formate  '1  ',  5</),  tSO,  'dimensions  of  errb  snkment  1  oad  ' ,  // ) 
63     formatetSO,  'embankment  hight  =  ',  f4.  1,  '  feet.  ',  /, 

1  t50,  'sideslope  =  ',f4.  1,  '  degrees  ',/, 

2  tSO,  'crest  width  =  ',  f5.  1,  '  feet.  ',  /, 

3  tSO,  'density  of  the  applied  load  =  ',  f5.  1,  '  pcf.  ') 

72     format(3(/), tl8, 'settlement  will  be  calculated  at  ',i2,     ^  ^ 

1  '  equally  spaced  positions  betueen  x=  ',f5. 1,  '  and  x-   , 

2  f5.  1,  '  feet.  '.  3(/)  ) 

74  format  (t45,  'layer  #',  15x,  'no.  of  d  i  vi  s  i  ons ',  // ) 

75  f ormat(t47,  i2,  t74,  i2) 

35     format  ('l',t25,  'x   (  feet )',  t40,  '  lay  er ',  t50,  '  strata  , 

1  t66,  'depth   (feet)  ',.t84,  'thickness  (feet)', 

2  tl05,  'strain  (  perc ent )  ' , 3 ( / )  ) 

210  formate  '1  ',  tSO,  'settlement  (  feet )',////- t55,  '  x  (feet)'.//) 

90  format  e  t26,  f6.  1  ) 

?1  format  (  '+ ' ,  t41 ,  i2,  t8a,  f  8.  1  ) 

92  format  e  '+ ' ,  t52,  i2,  t63,  f  7.  1  ) 

99  format  (  '+ ' ,  tl  10,  f  7.  2) 
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250    f  ormat  ( t25,  11  (f5.  1,  5x  )  ) 

305    format (//.  til,  'layer  '  ) 

310    format  (t  12,  12,  t25,  IKfS.  3,  5x  )  ) 

410    formate//,  til,  'total',  t25,  ll(f5.  3,  5x  )  ) 

end 

subroutine  stress  < x,  z,  hembnk,  slopel < wi d th ,  densty, s i gmaz ,  p i ) 
c 
c 

p  =  densty  *  hembnk 
c 

c  geometry 

c 
c 

slope    =    slopel    »    pi/ISO.  0 

a    =    hembnk/tan ( si  ope ) 

b    =    a    +    iuidth/2.  0 

xl    =    b    -    X 
c 

c  contribution  to  vertical  stress  from  right  half  of  the  load 

c 

r2  =  sqrt ( (width/2.  0  +  a  -  xl)*»2  +  2»*2) 

rl  =  sqrt<(xl  -  a)**2  +  2»«2) 

rO  =  sqrt(xl»*2  +  z*«2) 

betar  =  acos<<rl»«2  +  r2**2  -  (lui  d  th/2.  0)-»-«2)  /  (2,  0*r  l*r2)  ) 

gammar  =  aco5<<rl»*2  +  r0**2  -  a«*2) / ( 2, 0*r l*rO ) ) 

ciflr  =  betar  +  xl/a  »  gammar  -  2*( x 1-b ) /r2**2 
c 
c 

c  contribution  to  vert,  stress  from  the  left  side. 

c 
c 

rl  =  sqrt((u;idth  +  a  -  xl)**2  +  z»*2) 

x2  =  2.  0  »  a  +  width  -  xl 

rO  =  sqrt(x2K-*2  +  z-«-»2) 

betal  =  acos<<rl»«2  +r2**2  -  (ui d tb /2.  0 )»*2) / (2.  0*r 1 *r2 ) ) 

gammal  =  aco5((rl*«2  +  r0»*2  -  a**2  )  /  (  2  0*r  l-«-r0  )  ) 

cifll  =  betal  +  x2*gammal/a  -  2*( x2-b ) /r2#*2 
c 

sigmaz  =  p/pi*<ciflr  +  cifll) 

return 

end 
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APPENDIX  H 
TIME-RATE  OF  CONSOLIDATION  SETTLEMENT 

The  FORTRAN  IV  program  contained  hereafter  solves  the 
ujell  knou/n  Terzaghi  one-dimensional  consolidation  equation 


=vt:S=" 


by  means  of  the  finite-difference  approximation  discussed  in 
Chapter  IV.   Up  to  10  contiguous  layers  with  differing 
values  of  >z  _  _    may  be  input.   Any  distribution  of  initial 
excess  pore  pressure  may  be  input  because  the  solution  is 
obtained  numerically.   The  distribution  of  initial  excess 
pore  pressures  is  input  at  discrete  points.   These  points 
are  assumed  to  be  evenly  spaced  uithin  each  separate  layer. 
The  spacing  in  each  layer  is  determined  u;ith  the  following 
equation: 


ii:here 


^.|^....i'^^ 

Az.    =    — — — ''H 

1        cc . 


a  <  1/2 


i  =  layer  number 


th 


c  .  =  coef  f  ic  ierit  of  consol  idat  iori  of  thi='  i    lager 
'■-'  1 

t  =  time  increment  to  be  used  in  the  analysis 
liz .  =  spacing  of  nodes  in  the  i  '   layer. 


Whenever  possible  a  should  be  approximately  1/6. 
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Generally  the  solution  will  be  more  accurate  if  a  large 
number  of  nodal  points  is  chosen  in  each  layer. 
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User  Manual  -  Time  Rate  of  Settlement  Program 

All  input  is  unformatted. 

Units  are  feet,  pounds,  and  days  unless  otherwise 

spec  if  ied. 

1.  Read  in  the  number  of  contiguous  consolidating 
layers  (integer) 

2.  For  each  layer  read  in: 

a)  the  layer  number  (integer) 

b>  c_  of  the  layer  < ft  /day) 

c)  thickness  of  the  layer   (ft) 

d)  nodal  spacing  in  the  layer(ft) 

3.  Read  in  the  time  increment  at  which  excess  pore 
pressures  are  to  be  calculated  (days) 

^.   Read  in  the  values  of  the  initial  excess  pore 
pressures  at  each  of  the  nodes  from  top  to 
bottom  ( psf ) 

5.   Read  in  the  output  scheme:    (integer) 

a)  If  the  values  of  the  excress  pore 
pressure  at  the  nodes  are    desired, 
input  JL-   J^ue  to  format  restric- 
tions this  ouput  scheme  may  be  used 
only  if  15  or  fewer  nodes  are    used. 

b)  If  the  value  of  the  percent  consolidation 
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is  desired  in  each  layer,  input  2 

c)   If  the  value  of  the  percent  consolida- 
tion is  desired  for  a  one  layer  problem, 
input  0 

6.  Read  in  the: 

a>   code  for  the  drainage  condition  of 
the  top  boundary  (integer) 

b)  code  for  the  drainage  condition  of 
the  bottom  boundary  (integer) 

c)  depth  of  the  top  boundary  beneath 
the  ground  surface  (ft) 

If  a  boundary  is  drained,  set  its  drainage 
condition  code  =  0 

If  a  boundary  is  undrained,  set  its  drainage 
condition  code  =  1 

7.  Read  in  the  termination  limits: 

a)  the  maximum  number  of  years  at  which  the 
program  will  cease  to  calculate  the  ex- 
cess pore  pressures. 

b)  the  maximum  overall  percent  consolidation 
at  which  the  program  should  cease  to  calcu- 
late excess  pore  pressures. 

The  program  listing  is  contained  in  the  following 
pages. 
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INTEGER  RTOP»RBOTTMfCYCLESf I0UT»NSTRAT(10) 

REAL  U(101)fUINIT<101)»CV(10)»THICK(10) .DELZ(IO) r ALPHA ( 10) » 

1      Zll(20)fAINTL(10) f APRSNL(IO) »C0NS0L(10) 
C 
C 

C  PRINT  TITLE 

C 

WRITE  <6»1) 
WRITE  (6»2) 
C 
C  DIMENSIONS  OF  UNITS  IN  INPUT 

C  LENGTH  --  FEET 

C  PRESSURE  --  PSF 

C  TIME  —  DAYS 

C  EXCEPT  WHERE  SPECIFIED  OTHERWISE  AS  IN  THE  CASE 

C  OF  'NYEARS'  WHICH  IS  INPUT  IN  YEARS. 

C 

C  READ  IN  THE  NUMBER  OF  SEQUENTIAL  COMPRESSIBLE  LAYERS 

C 

READ  (5f*)  LAYERS 
C 

C  FOR  EACH  LAYERf  READ  INt 

C  1)  THE  LAYER  NUMBER  'II' 

C  2)  THE  LAYERS  COEFFICENT  OF  CONSOLIDATION  'CV(I1)'» 

C  3)  THICKNESS  'THICK(Il)' 

C  4)  THICKNESS  INCREMENT  'THICK  (ID' 

C 

DO  10  I  =  If  LAYERS 

READ  (5»*)  Il,CV<Il)fTHICK(Il>»DELZ(Il) 
10     CONTINUE 
C 

C  CALCULATE  THE  NUMBER  OF  STRATA  IN  EACH  LAYER  FOR 

C  COMPUTATIONAL  PURPOSES.   THE  USER  SHOULD  CHECK  THE  OUTPUT 

C  TO  INSURE  THAT  THE  PROGRAM  DIVIDES  EACH  LAYER  INTO  THE 

C  SAME  *  OF  STRATA  THAT  WAS  ASSUMED. 

C 

DO  20  I  =  1, LAYERS 
VALUEl  =  THICK(I)/DELZ(I> 
NSTRAT(I)  =  INT(VALUEl) 
VALUE2  =  FLOAT(NSTRAT(I>) 

IF  (VALUE2.lt. VALUED  NSTRAT(I)  =  NSTRAT(I)  +  1 
20     CONTINUE 
C 

C  READ  IN  THE  TIME  INCREMENT  'DELT'  TO  BE  USED  IN  THE  ANALYSIS 

C 

READ  (5f*)  DELT 
C 
C  CALCULATE  THE  VALUE  OF  'ALPHA'  FOR  EACH  LAYER 
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DO  30  I  =  If  LAYERS 
ALPHA(I)  =  CV(I)*DELT/DELZ<I)**2 
30     CONTINUE 
C 

C  CALCULATE  THE  TOTAL  *  OF  STRATA  'NS'  IN  ALL  OF  THE  LAYERS» 

C  AND  THE  NUMBER  OF  POINTS  'NPTS'  FOR  WHICH  THE  FINITE  DIFF- 

C  ERENCE  PROCEDURE  WILL  BE  PERFORMED. 

C 

NS  =  0 

DO  40  I  =  IrLAYERS 
NS  =  NS  +  NSTRAT(I) 
40     CONTINUE 

NPTS  =  NS  +  1 
C 

C  READ  IN  THE  VALUES  OF  THE  INITIAL  EXCESS  PORE  PRESSURE  WHICH 

C  ARE  ASSUhED  TO  BE  EQUAL  TO  THE  TOTAL  STRESS  CHANGE  AT  THE 

C  POINT  FOR  EACH  OF  THE  POINTS.   DO  NOT  REPEAT  VALUES  AT 

C  INTERLAYER  BOUNDARIES. 

C 
C 

READ  (5f*)  (UINIT<I).I=lfNPTS) 
C 
C 

C  PRINT  HEADINGS 

C 

WRITE  (6»55) 
C 

C  PRINT  INPUT  DATA 

C 

DO  70  I  =  If  LAYERS 

WRITE  (6»60)  IfTHICK(I)»NSTRAT(I)FDELZ(I) »CV(I) fALPHAd 
70     CONTINUE 
C 

C  CHECK  IF   0.0  <   ALPHA   <  0.5 

C 

DO  80  I  =  IfLAYERS 

IF  (ALPHA(I) .LE.0.5)  GO  TO  80 
WRITE  (6»75)  I 
STOP 
80     CONTINUE 
C 

C  READ  IN  OUTPUT  SCHEME 

C  OUTPUT  ONLY  THE  TIME  AND  THE  PERCENT  CONSOLIDATION 

C  SET   lOUT  =  0 

C  OUTPUT  THE  TIME,  THE  EXCESS  PORE  PRESSURE  AT  EACH  POINT 

C  AND  THE  PERCENT  CONSOLIDATION 

C  SET   lOUT  =  1 
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C 
C 

c 
c 
c 
c 
c 
c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c 
c 
c 


81 

82 

C 

C 

C 

C 

c 


95 
89 


OUTPUT  THE  TIME  .THE  PERCENT  CONSOLIDATION  IN  EACH  LAYER 
AND  THE  OVERALL  PERCENT  CONSOLIDATION 
SET  lOUT  =  2 

THE  USER  MUST  BE  CAREFUL  TO  NOT  USE  MORE  THAN  15  POIMTS  IN 
ORDER  TO  PREVENT  THE  OUTPUT  FORMATS  FROM  BLOWING  UP  WHEN 
USING  lOUT  =1 

READ  (5»*)  lOUT 

INPUT  BOUNDARY  CONDITIONS  OF  THE  TOP  AND  BOTTOM  OF  THE 
SOIL  LAYER  WHOSE  TIME  RATE  OF  SETTLEMENT  IS  BEING 
STUDIED. 

IF  THE  TOP  OR  BOTTOM  BOUNDARY  IS  DRAINED  THEN  INPUT 

RTOP  =  0   OR  RBOTTM  =  0   RESPECTIVELY. 
IF  THE  TOP  OR  BOTTOM  BOUNDARY  IS  UNDRAINEDf  THEN  INPUT 

RTOP  =  1  OR  RBOTTM  =  1   RESPECTIVELY. 

READ  (5,*)  RTOPfRBOTTMfZTOP 
INITIALIZATIONS 


CNS  =0.0 
T    =  0.0 

AINIT   =  0.0 

APRSNT  =0.0 

IF  (LAYERS 

DO  81 


8: 


,E0.1 )  GO  TO 
I=1»LAYERS 

AINTL(I)  =  0. 
APRSNL(I)  =  0. 
CONTINUE 
CONTINUE 

OUTPUT  THE  INITIAL  PORE  PRESSURES  AT  EACH  POINT  FROM  TOP 
TO  BOTTOM  AS  IT  WAS  INPUT.   THIS  FORMAT  REPEATS  THE 
PRESSURE  AT  INTER-STRATA  BOUNDARIES. 

IF  (lOUT.GT.O)  GO  TO  89 

NSTART  =  1 

WRITE  (6,85) 

DO  95  I  =  If  LAYERS 

NEND  =  NSTRAT(I)  +  NSTART 

WRITE  (6f90)  If <UINIT(J),J=NSTART»NEND) 

NSTART  =  NEND 

CONTINUE 

GO  TO  101 

IF  (lOUT.GT.l)  GO  TO  239 

WRITE  (6fl40) 

Zll(l)  =  ZTOP 

NP  =  1 


96 

CON 

97 

CONTINUE 

WRITE  (6»98) 

WRITE  (6f99) 

WRITE  (6fl50) 

WRITE  <6fl65) 

GO  TO  101 
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DO  97  I  =  1»LAYERS 

NENEi  =  NSTRAT(I) 
DO  96  J  =  1»NEND 
NP  =  NP  +  1 

Zll(NP)  =  ZlKNP-l)  +  DELZ(I) 
TINUE 

(Zll(NP)fNP=l»NPTS) 

Tf (UINIT(J)f J=1»NPTS) 
CNS 
GO  TO  101 

C 

c 

C  CALCULATE  THE  INITIAL  AREA  OF  THE  ISOCHRONE » I . E . » THE 

C  INTEGRAL  OF  THE  ORDINATES  OF  EXCESS  PORE  PRESSURE  * 

C  THE  DEPTH  INTERVAL  FOR  WHICH  THAT  VALUE  OF  PRESSURE 

C  IS  AN  AVERAGE  VALUE. 

C 

c 

239    WRITE  (6»339)  (  I » 1  =  1 , LAYERS ) 

101  NSTART  =  1 

DO  110  I  =  1»LAYERS 

NEND  =  NSTRAT(I)  +  NSTART 

NENDl  =  NEND  -  1 

DO  100  J  =  NSTARTfNEND 

IF  ( J. NE. NSTART)  60  TO  102 

AINIT  =  AINIT  +  UINIT(J)*DELZ< I)/2.0 

IF    (LAYERS. GE. 2)    AINTL(I)    =    AINTL(I)    +    UINIT ( J ) *DELZ ( I ) /2 . ( 

GlO    to    100 

102  IF  (J. GT. NENDl)  GO  TO  103 
AINIT  =  AINIT  +  UINIT( J)*DELZ(I) 

IF  (LAYERS. GE. 2)  AINTL(I)  =  AINTL(I)  +  UINIT ( J ) *DELZ ( I ) 
GO  TO  100 

103  AINIT  =  AINIT  +  UINIT ( J ) *DELZ ( I ) /2 . 0 

IF  (LAYERS. GE. 2)  AINTL(I)  =  AINTL(I)  +  UINIT ( J ) *DELZ ( I ) /2 . 
100    CONTINUE 

NSTART  =  NEND 

110  CONTINUE 
C 

C 

C  CHANGE  THE  EXCESS  PORE  PRESSURE  AT  UNCONFINED  BOUNDARIES 

C  TO  THE  AMBIENT  VALUE  =  0.5  *  INITIAL  EXCESS  PRESSURE. 

C 

C  ^ 

IF  (RTOP.EQ.O.AND.RBOTTM.EQ.l)  GO  TO  111 

IF  (RTOP.EQ.l.AND.RBOTTM.EQ.O)  GO  TO  112 

IF  (RTOP.EQ.O.AND.RBOTTM.EQ.O)  GO  TO  113 

IF  (RTOP.EQ.l.AND.RBOTTM.EQ.l)  GO  TO  114 

111  UINIT(l)  =  UINIT(l)/2.0 

APRSNT  =  AINIT  -  UINIT ( 1 ) *DELZ ( 1  ) /2  .  0 
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IF  (LAYERS. GE. 2)  APRSNL(l)  =  AINTL(l)  -  UINIT ( 1 ) *DELZ ( 1 ) /2 . 0 
GO  TO  120 

112  UINIT(NPTS)  =  UINIT(NPTS)/2.0 

APRSNT  =  AINIT  -  UINIT ( NPTS ) *DELZ ( LAYERS ) /2 . 0 

IF  (LAYERS. GE. 2)  APRSNL ( LAYERS)  =  APRSNL ( LAYERS )  -  UINIT(NPTS) 
1  *DELZ(LAYERS)/2.0 

GO  TO  120 

113  UINIT(l)  =  UINIT(l)/2.0 
UINIT(NPTS)  =  UINIT(NPTS)/2.0 

APRSNT  =  AINIT  -  ( UINIT ( 1 ) *DELZ ( 1 )  +  UINIT ( NPTS ) *DELZ ( LAYERS )) /2 . ( 

IF  (LAYERS. GE. 2)  APRSNL(l)  =  AINTL(l)  -  UINIT  ( 1)  *riELZ  ( 1 ) /2  .  0 

IF  (LAYERS. GE. 2)  APRSNL ( LAYERS )  =  APRSNL ( LAYERS )  -  UINIT(NPTS) 

GO  TO  120 
116    WRITE  (6rll7) 

STOP 
120    CNS  =  (1.0  -  APRSNT/AINIT)*100.0 

IF  (LAYERS. EG). 1)  GO  TO  126 

DO  125  I  =  If  LAYERS 

CONSOL(I)  =  (1.0  -  APRSNL(I)/AINTL(I))*100, 

125  CONTINUE 

126  CONTINUE 
C 

IF  (lOUT.EQ.O)  WRITE  (6fl30) 
IF  (lOUT.NE.O)  GO  TO  114 
WRITE  (6»145)  T»CNS 
GO  TO  115 

114  IF  (lOUT.NE.l)  GO  TO  115 

WRITE  (6»150)  T» (UINIT(J) ,J=1»NPTS) 

WRITE  (6»165)  CNS 
C 

C  INPUT  THE  ITERATION  LIMITS 

C  NYEARS  =  ♦  OF  YEARS  THE  PROGRAM  WILL  CALCULATE  EXCESS  PORE 

C  PRESSURES  AND  CONSOLIDATIONS. 

C  CMAX    =  THE  MAXIMUM  PERCENT  OF  CONSOLIDATION  THE  PROGRAM 

C  WILL  CALCULATE  PRESSURES  AND  CONSOLIDATIONS. 

C 

115  READ(5»*)  NYEARS»CMAX 
C 

CYCLES  =  NYEARS*365/INT(DELT) 
C 
C 
C 

C  FINITE  DIFFERENCE  LOOP 

C 
C 

C        CALCULATE  EXCESS  PORE  PRESSURES  AND  PERCENT  CONSOLIDATION 
C        FOR  EACH  TIME  STEP 
C 
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DO  200  ICYCLE  =  1»CYCLES 
T  =  T  +  DELT 
NSTART  =  1 
C 

C  CALCULATE  THE  EXCESS  PORE  PRESSURE  AT  ALL  THE  POINTS 

C  FOR  THE  ICYCLE'TH  CYCLE 

C 

DO  160  I  =  If  LAYERS 
NEND  =  NSTRAT(I)  +  NSTART 
C 

C  INTERNAL  POINTS  IN  EACH  LAYER 

C 

DO  155  J  =  NSTART » NEND 

IF  < (I.EQ.l.AND.J.EQ.l) )  GO  TO  155 

IF  (J.LT.NEND)  U(J)  =  ALPHA ( I ) * ( UINIT ( J+1 ) +UINIT < J-1 ) ) 
1  +  (1.0  -  2.0*ALPHA(I))*UINIT(J) 

155    CONTINUE 

NSTART  =  NEND 
160    CONTINUE 
C 

C  INTER  LAYER  BOUNDARY  POINTS 

C  IMPOSE  CONTINUITY  OF  FLOW  ACROSS  THE  BOUNDARY 

C 

NFACE  =  1 

NLl  =  LAYERS  -  1 

IF  (NLl .LE.O)  GO  TO  181 

DO  180  I  =  IfNLl 

NFACE  =  NFACE  +  NSTRAT(I) 

U(NFACE)  =  U(NFACE  +  1)  -  (U(NFACE  +  1)  -  U(NFACE  -  1) 
1  '  /(l.O  +  (CV(I+1)/CV(I) )»(DELZ(I)/DELZ(I+1) ) 

180  CONTINUE 
C 

C  CALCULATE  THE  PORE  PRESSURE  AT  THE  TOP  AND  BOTTOM  OF  THE 

C  SEQUENCE  OF  COMPRESSIBLE  LAYERS.  IT  IS  ASSUMED  THAT 

C  UNDRAINED  BOUNDARIES  MAY  BE  SIMULATED  WITH  A  MIRROR 

C  IMAGE  OF  THE  EXISTING  PORE  PRESSURES  ON  THE  OTHER  SIDE 

C  OF  THE  BOUNDARY. 

C 

C 

181  IF  (RTOP.EQ.O)  U(1)=0.0 

IF  (RTOP.EQ.l)  U(l)  =  2.0*ALPHA(1)*UINIT(2)  + 
1  (1.0  -  2.0*ALPHA( 1) )*UINIT(1) 

C 

IF  (RBOTTM.EQ.O)  U(NPTS)  =  0.0 

IF  (RBOTTM.EQ.l)   U(NPTS)  =  2 . 0*ALPHA ( LAYERS ) «UINIT ( NPTS-1 >  + 
1  (1.0  -  2.0*ALPHA(LAYERS) )*UINIT(NPTS) 
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C 

C  CALCULATE  THE  PERCENT  CONSOLIDATION 

C 

APRSNT  =0.0 

IF  (LAYERS. EQ.l)  60  TO  186 

DO  185  I  «  IfLAYERS 

APRSNL<I)  =  0. 

185  CONTINUE 

186  CONTINUE 

NSTART  =  1 

C 

C  CALCULATE  THE  PRESENT  ISOCHRONE 

C 

DO  195  I  =  1» LAYERS 
NEND  =  NSTRAT(I)  +  NSTART 
NENDl  =  NEND  -  1 
C 

DO  190  J  =  NSTARTfNEND 

IF  (J. NE. NSTART)  GO  TO  1181 
APRSNT  =  APRSNT  +  U ( J ) *DELZ ( I ) /2 . 0 
IF  (LAYERS. EQ. 2)  APRSNL(I)  =  APRSNL(I)  +  U(J>« 
1  DELZ(I)/2.0 

GO  TO  190 
1181  IF  (J. GT. NENDl)  GO  TO  182 

APRSNT  =  APRSNT  +  U(J>*DELZ(I) 

IF  (LAYERS. EQ. 2)  APRSNL(I)  =  APRSNL(I)  +  U(J)* 
1  DELZ(I) 

GO  TO  190 
192  IF  (J.EQ.NEND)  APRSNT  =  APRSNT  +  U( J ) *DELZ ( I ) /2 . 0 

IF  (J. EQ. NEND. AND. LAYERS. EQ. 2) 
1  APRSNL(I)=  APRSNL(I)  +  U ( J ) *DELZ ( I ) /2 . 0 

190    CONTINUE 

NSTART  =  NEND 
195    CONTINUE 

CNS  =  (1.0  -  APRSNT/AINIT)*100.0 
IF  (LAYERS. EQ.l)  GO  TO  1196 
DO  1195  I  =  1»LAYERS 

CONSOL(I)  =  (1.0  -APRSNL(I)/AINTL(I) )*100. 

1195  CONTINUE 

1196  CONTINUE 
C 

C  OUTPUT  THE  TIME»THE  PERCENT  CONSOLIDATION f 

C  AND  THE  EXCESS  PORE  PRESSURE  AT  EACH  NODE 

C 

IF  (I0UT.6T.0)  60  TO  196 

WRITE  (6»145)  T»CNS 

60  TO  197 
196    IF  (I0UT.6T.1)  60  TO  1197 

WRITE  (6»150)  T»(U(J) »J=1»NPTS) 

WRITE  (6»165)  CNS 

60  TO  197 
1197   WRITE  (6fl200)  T »( CONSOL ( I )> 1=1 r LAYERS ) 
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C 

c 
c 

197 
C 

c 
c 
c 


199 

200 

C 

300 

C 

C 

C 

1 


60 
75 


CHECK  TERMINATION  LIMITS 


IF 
IF 


(CNS.GT.CMAX)  GO  TO  300 
(ICYCLE.EQ. CYCLES)  GO  TO 


300 


RESET  THE  EXCESS  PORE  PRESSURES  AT  ALL  THE  POINTS 
FOR  THE  NEXT  TIME  CYCLE 

DO  199  I  =  IfNPTS 
UINIT(I)  =  U(I) 
CONTINUE 
CONTINUE 

CONTINUE 

FORMATS 


FORMAT 

FORMAT 
FORMAT 


FORMAT 


OF  SETTLEMENT  PROGRAM' »/» 
GOODMAN' f///) 


(1H1,30(/) »T50» 'TIME  RATE 
T52» 'PROGRAMMED  BY  MARTY 
(1H1»20(/) ) 

(T22»'LAYER  '» T31 .' THICKNESS '» T43 »' *  OF  STRATA'fT58» 
'STRATA  THICKNESS' »T78» 'COEFF.  OF  CONSOLIDATION ', T105» 
'ALPHA' >///) 
(T23»I2»T33fF5.1>T47f I2»T64rF5.2rT85f F8.5»T104»F6.4) 

LAYER  'rI2»'  HAS  A  VALUE  OF  ALPHA  MORE  THAN' 
'  EITHER  REDUCE  THE  TIME  INCREMENT  DELT',/,T20 


FORMAT  ( ' 
'  0.5' f /, 
'OR' »/» 
'  REDUCE 
'  VALUES 


THE  *  OF  STRATA  IN  THE  LAYERS  WITH  CALCULATED' 
OF  ALPHA  GREATER  THAN  0.5' »/f 
5  '  THIS  MAY  BE  ACCOMPLISHED  BY  INCREAIN6  DELZ') 

85     FORMAT  (IHl , T5 »' LAYER '» T20. ' INITIAL  EXCESS  PORE  PRESSURES'.//) 

90     FORMAT  (  5X . 12 . 5X . 18 ( F5 . 0 f IX ) ) 

98  FORMAT  (lOX . 1 8 ( IX f F6 . 0 ) ) 

99  FORMAT  (2(/)) 

117    FORMAT  ('  THE  TOP  OF  AND  BOTTOM  BOUNDARIES  ARE  UNDRAINED ' . / » 

1  '  NO  COSOLIDATION  WILL  OCCUER') 

130    FORMAT  ( IHl . 5 ( / ) f T43 , ' TIME  ( DAYS )'. T67 »' PERCENT  CONSOLIDATION'/// 
140    F0RMAT(1H1»5(/) .T2» 'TIME  ( DAYS )'» T50 .' EXCESS  PORE  PRESSURE  (PSF)' 

1        T120. 'CONSOLIDATION' »/.T59. 'Z  (FEET)'./) 
145    FORMAT  ( T40 . FIO . 1 . 20X . FIO . 2 ) 
150    FORMAT  ( F8 . 1 . 2X . 15( IX . F6 . 0 ) ) 
165    FORMAT  ( 1H+ . Tl 19 . FIO . 2 ) 

339    FORMAT (1  HI.///, Til. 'TIME' .T50. ' LAYERS ',//. T25 . 10 ( 14 . 6X) .//) 
1200   FORMAT <T5. FIO, 1.8X. 10 (F8.1.2X)) 
END 


